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Chapter  6 

Multiple  Target  Tracking 


“Although  this  may  seem  a  paradox,  all  exact  science  is 
dominated  by  the  idea  of  approximation." 

Bertrand  Russell,  The  Scientific  Outlook,  1931 

Abstract  Multitarget  tracking  intensity  filters  are  closely  related  to  imaging 
problems,  especially  PET  imaging.  The  intensity  filter  is  obtained  by  three 
different  methods.  One  is  a  Bayesian  derivation  involving  target  prediction 
and  information  updating.  The  second  approach  is  a  simple,  compelling, 
and  insightful  intuitive  argument.  The  third  is  a  straightforward  application 
of  the  Shepp-Vardi  algorithm.  The  intensity  filter  is  developed  on  an  aug¬ 
mented  target  state  space.  The  PHD  filter  is  obtained  from  the  intensity  filter 
by  substituting  assumed  known  target  birth  and  measurement  clutter  inten¬ 
sities  for  the  intensity  filter's  predicted  target  birth  and  clutter  intensities. 

To  accommodate  heterogeneous  targets  and  sensor  measurement  models, 
a  parameterized  intensity  filter  is  developed  using  a  marked  PPP  Gaussian 
sum  model.  Particle  and  Gaussian  sum  implementations  of  intensity  fil¬ 
ters  are  reviewed.  Mean-shift  algorithms  are  discussed  as  a  way  to  extract 
target  state  estimates.  Grenander's  method  of  sieves  is  discussed  for  regu¬ 
larization  of  the  multitarget  intensity  filter  estimates.  Sources  of  error  in  the 
estimated  target  count  are  discussed.  Finally,  the  multisensor  intensity  filter 
is  developed  using  the  same  PPP  target  models  as  in  the  single  sensor  fil¬ 
ter.  Both  homogeneous  and  heterogeneous  multisensor  fields  are  discussed. 
Multisensor  intensity  filters  reduce  the  variance  of  estimated  target  count  by 
averaging. 


Multitarget  tracking  in  clutter  is  a  joint  detection  and  estimation  problem. 
It  comprises  two  important  inter-related  tasks.  One  initiates  and  terminates 
targets,  and  the  other  associates,  or  assigns,  data  to  specific  targets  and  to 
clutter.  The  MHT  (Multiple  Hypothesis  Tracking)  formulation  treats  both 
tasks.  It  is  well  matched  to  the  target  physics  and  to  the  sensor  signal  proces¬ 
sors  of  most  radar  and  sonar  systems.  Unfortunately,  exact  MHT  algorithms 
are  intractable  because  the  number  of  measurement  assignment  hypotheses 
grows  exponentially  with  the  number  of  measurements.  These  problems  are 


151 


5-2 


RTO-EN-SET-1 57(2010) 


Poisson  Point  Processes:  Imaging,  Tracking,  and  Sensing 


152  6  Applications:  Multiple  Target  Tracking 

aggravated  when  multiple  sensors  are  used.  Circumventing  the  computa¬ 
tional  difficulties  of  MHT  requires  approximation. 

Approximate  tracking  methods  based  on  PPP  models  are  the  topics  of 
this  chapter.  They  show  much  promise  in  difficult  problems  with  high  target 
and  clutter  densities.  The  key  insight  is  to  model  the  distribution  of  targets 
in  state  space  as  a  PPP,  and  then  use  a  filter  to  update  the  defining  parameter 
of  the  PPP  —  its  intensity.  To  update  the  intensity  is  to  update  the  PPP.  The 
intensity  function  of  the  PPP  approximation  characterizes  the  multiple  target 
tracking  model.  This  important  point  is  discussed  further  in  Section  6.1.1. 

The  PPP  intensity  model  uses  an  augmented  state  space,  <S+.  This  enables 
it  to  estimate  target  birth  and  measurement  clutter  processes  on-line  as  part 
of  the  filtering  algorithm. 

Three  approaches  to  the  intensity  filter  are  provided.  The  first  and  most 
rigorous  is  a  Bayesian  derivation  given  in  Appendix  D.  The  relationship 
between  this  approach  and  the  "first  moment  intensity"  of  the  posterior 
point  process  is  shown.  The  second  approach  is  a  short  but  extraordinarily 
insightful  derivation  that  is  ideal  for  readers  who  wish  to  avoid  the  Bayesian 
analysis,  at  least  on  a  first  reading.  The  third  approach  is  based  on  the 
connection  to  PET  and  the  Shepp-Vardi  algorithm.  The  PET  interpretation 
contributes  significantly  to  understanding  the  PPP  target  model.  A  special 
case  of  the  intensity  filter  is  the  well  known  PHD  filter.  It  is  obtained  by 
assuming  a  known  target  birth-death  process,  a  known  measurement  clutter 
process,  and  restricting  the  intensity  filter  to  the  non-augmented  target  state 
space  S. 

Implementation  issues  are  discussed  in  Section  6.3.  Current  approaches 
use  either  particle  or  Gaussian  sum  methods.  An  image  processing  method 
called  the  mean  shift  algorithm  is  well  suited  to  point  target  estimation, 
especially  for  particle  methods.  Observed  information  matrices  ( cf.  Section 
4.7)  are  proposed  as  surrogates  for  the  error  covariance  matrices  widely  used 
in  single  target  Bayesian  filters.  The  underlying  statistical  meaning  of  OIM 
estimates  is  as  yet  unresolved. 

Several  sources  of  error  in  the  estimated  target  count  are  discussed  in 
Section  6.4.  Target  count  estimates  are  sensitive  to  the  probability  of  target 
detection.  This  function,  P^(x),  depends  on  target  state  x  at  measurement 
time  4-  It  varies  over  time  because  of  slowly  changing  sensor  characteristics 
and  environmental  conditions.  Monitoring  these  and  other  factors  that  affect 
target  detection  probability  is  not  typically  considered  part  of  the  tracking 
problem. 

Several  areas  of  on-going  research  are  also  discussed.  One  is  a  Gaussian 
sum  PPP  filter  that  enables  heterogeneous  target  motion  and  sensor  mea¬ 
surement  models  to  be  used  in  an  intensity  filter  setting.  See  Section  6.2.2  for 
details.  Another  is  the  multiple  sensor  intensity  filter  described  in  Section  6.5. 
This  filter  is  relies  on  the  validity  of  the  PPP  approximation  for  every  sensor. 
It  reduces  the  variance  of  the  target  count  by  averaging  over  the  sensors. 
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so  that  the  variance  of  estimated  target  count  decreases  with  the  number  of 
sensors. 


6.1  Intensity  Filters 

6.1.1  PPP  Model  Interpretation 

The  points  of  a  realization  of  a  PPP  on  the  target  state  space  are  a  poor  rep¬ 
resentation  of  the  physical  reality  of  a  multiple  target  state.  This  is  especially 
easy  to  see  when  exactly  one  target  is  present,  for  then  ideally 

J  A(x)  dx  =  l.  (6.1) 

From  (2.4),  the  probability  that  a  realization  of  the  PPP  has  exactly  one  point 
target  is 

VN(n  =  l)=e~1  ~  37%. 

Hence,  63%  of  all  realizations  have  either  no  target  or  two  or  more  targets. 
Evidently,  realizations  of  the  PPP  seriously  mismodel  this  simple  tracking 
problem.  The  problem  worsens  with  increasing  the  target  count:  if  exactly  n 
targets  are  present,  then  the  probability  that  a  realization  has  exactly  n  points 
is  e~nnn  / n\  ~  ( 27in )-1/2  — >  0  as  n  — »  oo.  Evidently,  PPP  realizations  are  poor 
models  of  real  targets. 

One  interpretation  of  the  PPP  approximation  is  that  the  critical  element 
of  the  multitarget  model  is  the  intensity  function,  not  the  PPP  realizations. 
The  shift  of  perspective  means  that  the  integral  (6.1)  is  the  more  physically 
meaningful  quantity.  Said  another  way,  the  concept  of  expectation,  or  en¬ 
semble  average  over  realizations,  corresponds  more  closely  to  the  physical 
target  reality  than  do  the  realizations  themselves. 

A  huge  benefit  comes  from  accepting  the  PPP  approximation  to  the  mul¬ 
tiple  target  state  —  exponential  numbers  of  assignments  are  completely 
eliminated.  The  PPP  approximation  finesses  the  data  assignment  problem 
by  replacing  it  with  a  stochastic  imaging  problem,  and  the  imaging  problem 
is  easier  to  solve.  It  is  fortuitous  that  the  imaging  problem  is  mathematically 
the  same  problem  that  arises  in  PET;  see  Section  5.2.  The  "at  most  one  mea¬ 
surement  per  target"  rule  for  tracking  corresponds  in  PET  to  the  physics  — 
there  is  at  most  measurement  per  one  positron-electron  annihilation. 

Analogies  are  sometimes  misleading,  but  consider  this  one:  In  the  lan¬ 
guage  of  thermodynamics,  the  points  of  PPP  realizations  are  microstates. 
Microstates  obey  the  laws  of  physics,  but  are  not  directly  observable  without 
disturbing  the  system  state.  Physically  meaningful  quantities  (such  as  tem¬ 
perature,  etc.)  are  ensemble  averages  over  the  microstates.  In  the  PPP  target 
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model,  the  points  of  a  realization  are  thus  "microtargets",  and  microtargets 
obey  the  same  target  motion  and  measurement  models  as  real  targets.  The 
ensemble  average  over  the  PPP  microtargets  yields  the  target  intensity  func¬ 
tion.  The  language  of  microtargets  is  helpful  in  Section  6.5  on  multisensor 
intensity  filtering,  but  it  is  otherwise  eschewed  in  this  chapter. 


6.1.2  Predicted  Target  and  Measurement  Processes 

Formulation 

Standard  filtering  notation  is  adopted,  but  modified  to  accommodate  PPPs. 
The  general  Bayesian  filtering  problem  is  reviewed  in  Appendix  C.  Let  S  = 
R"'  denote  the  «v-dimensional  single  target  state  space.  The  augmented 
space  is  5+  =  SUcp,  where  (b  represents  a  "clutter  target"  state.  Clutter  targets 
account  for  data  not  assigned  to  real,  physical  targets.  The  augmented  space 
<S+  is  discussed  in  Section  2.12. 

The  single  target  transition  function  from  time  tk_i  to  time  tk,  denoted  by 
Wk_dy  I  A')  =  PEk\Ek_i  (y  I  x)r  is  assumed  known  for  all  x,yeS+.  The  augmented 
state  space  enables  both  target  initiation  and  termination  to  be  incorporated 
directly  into  as  specialized  kinds  of  state  transitions.  Novel  aspects  of 
the  transition  function  are: 

•  Wk_i((p\x)  is  the  probability  that  a  target  at  x  £  S  terminates; 

•  T'k-i (y  I  (p)  is  the  likelihood  that  a  new  target  initiates  at  y  e  S;  and 

•  W  ,  {(b  |  (p)  is  the  probability  that  a  target  in  <p  remains  in  <p. 

The  augmented  state  cp  is  also  used  to  account  for  measurement  clutter,  that 
is,  for  data  that  do  not  correspond  to  real  targets. 

The  multitarget  state  at  time  tk  is  Sk.  It  is  a  point  process  on  »S+,  but  it 
is  not  a  PPP  on  *S+.  Nonetheless,  Ek  is  approximated  by  a  PPP  to  "close 
the  loop"  after  the  Bayesian  information  update.  The  multitarget  state  is  a 
realization  E,k  of  Ek.  If  £,k  =  ( n ,  j.t] ,  ...,xn}),  then  every  point  X;  is  either  a 
point  in  S  or  is  cb.  It  is  stressed  that  repeated  occurrences  of  cp  are  allowed  in 
the  list  (xi,...,x„)  to  account  for  clutter. 

The  measurement  at  time  tk  is  Tk.  It  is  a  point  process  on  the  (nonaug- 
mented)  space  T  =  R"z,  where  nz  is  the  dimension  of  a  sensor  measurement. 
The  measurement  data  set 


vk  =  ( m ,  {zi,...,zm})e&{T) 

is  a  realization  of  Tk.  The  pdf  of  a  point  measurement  zeT  conditioned  on  a 
target  in  state  x  e  S+  at  time  tk  is  the  measurement  pdf  pk(z  \  x).  The  only  novel 
aspect  of  this  pdf  is  that  pk(z  \  (p)  is  the  pdf  that  z  is  a  clutter  measurement. 

The  Bayesian  posterior  multitarget  state  point  process  conditioned  on  the 
data  V\, ...,  vk~i  is  approximated  by  a  PPP.  Denote  this  PPP  by  Ek-i\k~i-  The 
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Intensity:  Jk  /x) 


Current  Multi-Target  PPP 


Predicted  Multi-Target  PPP 

1- 


PPP  Thinning 


Probability 
of  Detection 


New  Measurement  Set 


Pass 


Fail 


PET  Iteration 


Update 
Detected 
Target  PPP 

IZZ 


PPP  Superposition 


Add  intensities 

■ 


Updated  Multi-Target  PPP 


Updated  Intensity:  /^(x) 


Fig.  6.1  Block  diagram  of  the  Bayes  update  of  the  intensity  filter  on  the  augmented  target 
state  space  «S+.  One  step  of  the  Shepp-Vardi  algorithm  is  used  in  the  "PET  Iteration"  block. 
Because  the  null  state  tp  is  part  of  the  state  space,  target  birth  and  measurement  clutter 
estimates  are  intrinsic  to  the  predicted  target  and  predicted  measurement  steps.  The  same 
block  diagram  holds  for  the  PHD  filter  on  the  nonaugmented  space  S. 


intensity  of  Hfc-i|ic-i  is  fk-i\k-i(x)>  x  e  ^+.  Let  S^_ i  denote  the  predicted  PPP 
at  time  4-  Its  intensity  is  denoted  by  fk\k-i(x),  and  it  is  the  integral  of  the 
intensity  fk-i\k-i(x)  of  Sj_i|j._i,  as  seen  in  Section  2.11.2,  Eqn.  (2.86). 

The  goal  is  to  update  the  predicted  PPP,  Sj.ij._i,  with  the  measurement  data 
Vk-  The  information  updated  point  process  is  not  a  PPP,  so  it  is  approximated 
by  a  PPP.  Let  Sj.|j.  denote  the  approximating  PPP,  and  let  its  intensity  be 
fk\k(x)- 

Fig.  6.1  outlines  the  steps  of  the  intensity  filter.  The  discussion  below 
walks  through  the  steps  in  the  order  outlined. 
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Target  Motion  and  the  Bernoulli  Split 

The  first  of  these  steps  accounts  for  target  motion  and  predicts  the  intensity 
at  the  next  time  step.  The  input  is  a  PPP  with  intensity  fk-i\k-]  (x)>  so  the 
transition  procedure  yields  an  output  process  that  is  also  a  PPP,  as  is  seen  in 
Section  2.11.1.  Let  fk\k-i(x)  denote  the  intensity  of  the  output  PPP.  Adapting 
(2.83)  to  5+  gives 


fk\k-i(x)  =  f  ^{x\y)fk-i\k-i(y)dy,  (6.2) 

JS+ 

where  the  integral  over  <S+  is  defined  as  in  (2.97). 

Target  motion  is  followed  by  a  Bernoulli  thinning  procedure  using  the 
probability  of  the  sensor  detecting  a  target  at  the  point  x,  denoted  P^(x).  This 
probability  is  state  dependent  and  assumed  known.  The  input  PPP  intensity 
is  fk\k-i(x)-  As  seen  in  Section  2.9.2,  thinning  splits  it  into  two  PPPs  -  one 
for  detected  targets  and  the  other  for  undetected  targets,  denoted  by /^^(x) 

and  !  (x),  respectively.  These  PPPs  are  independent  (see  Section  2.9.2) 
and,  from  (2.56),  and  their  intensities  are 

fm-  ito  =  P?(x)fw c-iW 


and 


/^-iW  =  (1-p°W)/^-iW. 

Both  branches  in  Fig.  6.1  are  now  subjected  to  an  information  update. 


Predicted  Measurement  PPP,  and  Why  It  Is  Important 

As  seen  in  Section  2.11.2,  the  predicted  measurement  process  is  a  PPP.  Its 
intensity  is 


Ajt|fc_i(z)  =  pk(z\x)P,P(x)fk\k-i(x)dx,  for  zeT,  (6.3) 

JS+ 

as  is  seen  from  (2.86).  The  measurement  PPP  is  a  critical  component  of  the 
intensity  filter  because,  as  is  seen  in  Eqn.  (6.9),  it  weights  the  individual  terms 
in  the  sum  that  comprises  the  filter. 

Another  way  to  see  the  importance  of  A^-^z)  is  to  recall  the  classical 
single-target  Bayesian  tracking  problem.  The  standard  Bayesian  formulation 
gives 
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p(x\z)  = 


p{z\x)p(x) 

p(z) 


(6.4) 


where  the  denominator  is  a  scale  factor  that  makes  the  left  hand  side  a 
true  pdf.  It  is  very  easy  to  ignore  p(z)  in  practice  because  the  numerator 
is  obtained  by  multiplication  and  the  product  scaled  so  that  it  is  a  pdf. 
When  multiple  conditionally  independent  measurements  are  available,  the 
conditional  likelihood  is  a  product  and  it  is  the  same  story  again  for  the  scale 
factor.  However,  if  the  posterior  pdfs  are  slimmed,  not  multiplied,  the  scale  factor 
must  be  included  for  the  individual  terms  to  be  comparable.  Such  is  the  case 
with  the  intensity  filter:  the  PPP  model  justifies  adding  Bayesian  posteriors 
instead  of  multiplying  them,  and  the  scale  factors  are  crucial  to  making  the 
sum  meaningful. 

The  scale  factor  clearly  deserves  a  respectable  name,  and  it  has  one.  It  is 
called  the  partition  function  in  statistical  physics  and  the  machine  learning 
communities. 


6.1.3  Information  Updates 

It  is  seen  in  Appendix  C  that  the  mathematically  correct  information  update 
procedure  is  to  apply  the  Bayesian  method  to  both  the  detected  and  unde¬ 
tected  target  PPPs  to  evaluate  their  posterior  pdfs.  These  pdfs  are  defined  on 
the  event  space  <?(»S+).  If  the  posterior  pdfs  have  the  proper  form,  then  the 
posterior  point  processes  are  PPPs  and  are  characterized  by  their  intensity 
functions  on  S+. 

The  information  update  of  the  undetected  target  PPP  is  the  Bayesian 
updated  process  condition  on  no  target  detection.  The  posterior  point  process 
is  identical  to  the  predicted  target  point  process.  It  is  therefore  a  PPP  whose 
intensity,  denoted  by  f^k(x),  is  given  by 

/£(*)  =  fSc-  iM 

=  (l  -  P°(x))fk[k-f  x).  (6.5) 

This  brings  the  right  hand  branch  of  Fig.  6.1  to  the  superposition  stage,  i.e., 
the  block  that  says  "Add  Intensities". 

The  left  hand  branch  is  more  difficult  because,  as  it  turns  out,  the  infor¬ 
mation  updated  detected  target  point  process  is  not  a  PPP.  This  is  a  serious 
dilemma  since  it  is  highly  desirable  both  theoretically  and  computationally 
for  the  filter  recursion  to  remain  a  closed  loop.  The  posterior  point  process 
of  the  detected  targets  is  therefore  approximated  by  a  PPP.  Three  methods 
are  given  for  obtaining  this  approximation. 

The  first  method  is  a  Bayesian  derivation  of  the  posterior  density  of  the 
point  process  on  the  event  space  S’(S+),  followed  by  a  "mean  field"  ap- 
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proximation.  Details  about  Bayesian  filters  in  general  are  given  in  Appendix 
C.  The  Bayesian  derivation  is  mathematically  rigorous,  but  not  particularly 
insightful. 

To  gain  intuition,  one  need  look  no  further  than  to  the  second  method. 
While  not  rigorous,  it  is  intuitively  appealing  and  very  convincing.  The  per¬ 
spective  is  further  enriched  by  the  third  method.  It  shows  a  direct  connection 
between  the  information  update  of  the  Bayesian  filter  and  the  Shepp-Vardi 
algorithm  for  PET  imaging.  This  third  method  also  poses  an  interesting 
question  about  iterative  updates  of  a  Bayesian  posterior  density. 


First  Method:  Bayesian  Derivation  and  Mean  Field  Approximation 

The  pdf  of  the  Bayesian  posterior  point  process  for  detected  targets  is  defined 
on  the  event  space  <%’(S+).  The  Bayesian  posterior,  or  information  updated, 
pdf  is  defined  on  this  complex  event  space.  The  derivation  is  straightforward 
and  a  delight  for  Bayesians.  It  is  relatively  long  and  interferes  with  the  flow 
of  the  discussion,  so  it  is  given  in  Appendix  D  where  it  can  be  read  at  leisure. 
The  main  points  are  outlined  here,  so  readers  suffer  little  or  no  loss  of  insight 
by  skipping  it  on  a  first  reading.  Specific  equation  references  are  provided 
here  so  the  precision  of  the  mathematics  is  not  lost  in  the  flow  of  words. 

The  Bayesian  derivation  explicitly  incorporates  the  "at  most  one  measure¬ 
ment  per  target"rule  into  the  measurement  likelihood  function.  It  imposes 
this  constraint  via  the  measurement  conditional  pdf  ( cf.  Eqn.  (D.5)).  This  pdf 
sums  over  all  possible  assignments  of  the  given  data  to  targets.  Because  of  the 
augmented  space  S+ ,  a  clutter  measurement  is  accounted  for  by  assigning  it 
to  a  target  with  state  <ft. 

The  usual  Bayesian  update  (cf.  Eqn.  (D.6))  leads  to  the  pdf  of  the  Bayesian 
posterior  point  process  (cf.  Eqn.  (D.10)).  This  pdf  uses  the  facts  that  the  a  priori 
target  process  is  a  PPP  with  intensity  f^k_^(x)  and  that  the  predicted  measure¬ 
ment  process  is  a  PPP  with  intensity  Ak\k-\  (z)  given  by  (6.3).  The  posterior 
pdf  is  computationally  intractable  except  in  reasonably  small  problems.  In 
any  event,  inspection  of  the  posterior  pdf  clearly  reveals  that  it  does  not  have 
the  form  of  a  PPP  pdf.  Approximating  the  posterior  point  process  with  a  PPP 
is  the  next  step. 

The  Bayesian  posterior  pdf  is  replaced  by  a  mean  field  approximation,  a 
widely  used  method  of  approximation  in  machine  learning  and  statistical 
physics  problems.  This  approximation  is  the  product  of  the  one  dimensional 
marginal  pdfs  of  the  posterior  pdf.  The  marginal  pdfs  are  identical.  The 
intensity  function  of  the  approximating  PPP  is  therefore  taken  to  be  propor¬ 
tional  to  the  marginal  pdf  (cf.  Eqn.  (D.13)).  The  appropriate  scale  factor  is 
a  constant  determined  by  maximum  likelihood.  In  essence,  the  mean  field 
approximation  (cf.  Eqn.  (D.14))  is  proportional  to  the  intensity  function  of 
the  approximating  PPP. 
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The  posterior  detected  target  point  process  is  a  special  case  of  the  more 
general  class  of  finite  point  processes.  The  general  theory  of  these  processes 
dates  to  the  1950s.  An  excellent  text  on  the  topic  is  [16,  Chap.  5].  Finite  point 
process  theory  leads  via  the  so-called  Janossy  densities  to  the  first  moment 
intensity  function  approximation.  Janossy  densities  are  merely  the  joint  pdfs 
of  the  points  of  the  process  (conditioned  on  their  number).  The  first  moment 
is  a  sum  over  all  Janossy  densities.  In  the  tracking  application,  only  one 
Janossy  density  is  nonzero,  and  it  is  the  likelihood  function  of  updated  point 
process.  Consequently,  and  happily,  the  first  moment  is  given  by  a  sum 
comprising  only  one  term.  The  mean  field  approximation  is  very  closely 
related  to  the  first  moment  intensity  function.  Further  details  are  provided 
in  the  last  section  of  Appendix  D. 


Second  Method:  Expected  Target  Count 


Let  pk(  ■  |  x),  x  e  S+,  denote  the  conditional  pdf  of  a  measurement  in  the  mea¬ 
surement  space  T ,  so  that 


J'  pk(z  |  x)  dz  =  1 ,  for  all  x  e  S+ . 


(6.6) 


The  special  case  pk(z  |  cp)  is  the  pdf  of  a  data  point  z  conditioned  on  the  target 
state  (p ,  that  is,  the  pdf  of  z  given  that  it  is  clutter.  The  predicted  intensity  at 
time  tk  of  the  target  point  process  is  a  PPP  with  intensity  P®(x)fk \k-\(x).  The 
intensity  f^k(x)  is  the  intensity  of  a  PPP  that  approximates  the  information 
updated,  or  Bayes  posterior,  detected  target  process. 

The  measured  data  at  time  tk  are  mk  points  in  a  measurement  space  T . 
Denote  these  data  points  by  Zk  =  (zj, .  ..,zntk). 

The  information  update  of  the  detected  target  PPP  is  obtained  intuitively 
as  follows.  The  best  current  estimate  of  the  probability  that  the  point  mea¬ 
surement  Zj  originated  from  a  physical  target  with  state  x  e  S  in  the  infinites¬ 
imal  dx  is 


Pk(zj\x)Pk(x)fk\k-i(x)dx 

A-k\k-l(zj) 


(6.7) 


where  the  denominator  is  found  using  (6.3).  Similarly,  the  probability  that  z; 
originated  from  a  target  with  state  cp  is 


Pk(Zj\<p)Pk((p)fM-l(cp) 

hk\k-l(Zj) 


(6.8) 
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Because  of  the  "at  most  one  measurement  per  target"  rule,  the  sum  of  the 
ratios  over  all  measurements  Zj  is  the  estimated  number  of  targets  at  x,  or 
targets  in  cp,  that  generated  a  measurement. 

The  estimated  number  of  targets  at  x  £  S  is  set  equal  to  the  expected 
number  of  targets  conditioned  on  the  data  vk,  namely  fuk(x)  dx.  Cancelling 
dx  gives 


&  w  =  £ 

7=1 


Pk(z j  \x)?k  (x)  fk\k-l{x) 
hk\k-l(Zj) 


(6.9) 


Eqn.  (6.9)  holds  for  all  x  £  S+,  not  just  for  x  £  S. 

The  expected  target  count  method  makes  it  clear  that  the  expected  number 
of  detected  targets  in  any  given  set  *R  c  S  is  simply  the  integral  of  the  posterior 
pdf 


E  [Number  of  detected  targets 


(6.10) 


Similarly,  the  expected  number  of  targets  in  state  c p  is  the  posterior  intensity 
evaluated  at  cp,  namely,  fuk(<p)-  The  predicted  measurement  process  with 
intensity  A/t|jt-i(z)  is  a  vital  part  of  the  intensity  filter. 


Third  Method:  Shepp-Vardi  Iteration 

The  PET  model  is  interesting  here.  The  measurement  data  and  multiple 
target  models  are  interpreted  analogously  so  that: 

•  The  target  state  space  b>+  corresponds  to  the  space  in  which  the  radioiso¬ 
tope  is  absorbed. 

•  The  measured  data  Zk  c  T  correspond  to  the  measured  locations  of  the 
annihilation  events.  As  noted  in  the  derivation  of  Shepp-Vardi,  the  mea¬ 
surement  space  7~  need  not  be  the  same  as  the  state  space. 

•  The  posterior  target  intensity  f^k(x)  corresponds  to  the  annihilation  event 
intensity. 

The  analogy  makes  the  targets  mathematically  equivalent  to  the  distribution 
of  (hypothetical)  positron-electron  annihilation  events  in  the  state  space. 

Under  the  annihilation  event  interpretation,  the  information  update  (6.9) 
of  the  detected  target  process  is  given  by  the  Shepp-Vardi  algorithm  for  PET 
using  PPP  sample  data.  The  EM  derivation  needs  only  small  modifications 
to  accommodate  the  augmented  state  space.  Details  are  left  to  the  reader. 
The  /i-th  iteration  of  the  Shepp-Vardi  algorithm  is,  from  (5.18), 
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mk 

®-)M  =  /£ww  £ 


jri  Js.W(Zjls)/5(s)w<is' 


Pk(Zj\x) 


(6.11) 


where  the  predicted  intensity  fR.(x)*®  -  /^_1(x)  =  Pk(x)fk\k-i(x)  initializes 
the  algorithm.  The  first  iteration  of  this  version  of  the  Shepp-Vardi  algorithm 
is  clearly  identical  to  the  Bayesian  information  update  (6.9)  of  the  detected 
target  process.  The  second  and  higher  iterations  are  not  Bayesian  intensity 
estimates. 

The  Shepp-Vardi  iteration  converges  to  an  ML  estimate  of  the  target  state 
intensity  given  only  data  at  time  t k.  It  is  independent  of  the  data  at  times 
t\,  fjt_i  except  insofar  as  the  initialization  influences  the  ML  estimate.  In 
other  words,  the  iteration  leads  to  an  ML  estimate  of  an  intensity  that  does 
not  include  the  effect  of  a  Bayesian  prior.  The  problem  lies  not  in  the  PET 
interpretation  but  in  the  pdf  of  the  data.  To  see  this  it  suffices  to  observe 
that  the  parameters  of  the  pdf  (5.7)  are  not  constrained  by  a  Bayesian  prior 
and,  consequently,  the  Shepp-Vardi  algorithm  converges  to  an  estimate  that 
is  similarly  unconstrained.  It  is,  moreover,  not  obvious  how  to  impose  a 
Bayesian  prior  on  the  PET  parameters  that  does  not  disappear  in  the  small 
cell  limit. 


6.1.4  The  Final  Filter 

Superposing  the  PPP  approximation  of  the  detected  target  process  and  the 
undetected  target  PPP  gives 


fk\k(x)  =  f^k(x)  +  ffik(x),  «5+ 


as  the  updated  intensity  of  the  PPP  approximation  to  Ekik. 

The  intensity  filter  comprises  equations  (6.2),  (6.3),  and  (6.12).  The  first 
two  equations  are  more  insightful  when  written  in  traditional  notation.  Ex¬ 
panding  the  discrete-continuous  integral  (6.2)  gives 


(6.13) 


where  the  predicted  target  birth  intensity  is 

bk(x)  =  Wk-i(x\cp)  fk—i\k—i((p). 


(6.14) 


Also,  from  (6.3), 


5-12 


RTO-EN-SET-1 57(2010) 


Poisson  Point  Processes:  Imaging,  Tracking,  and  Sensing 


162  6  Applications:  Multiple  Target  Tracking 

Table  6.1  Intensity  Filter  on  the  State  Space  S+  =  S  U  p 
INPUTS: 

Data:  zj:m  =  ,zm)  <zT  at  time  4 
Probability  of  target  detection:  P°(x)  at  time  tk 

OUTPUT:  (fk\k(x),  fk\k(<p )l  =  IntensityFilter  [/£_!!*_!  (x),  fk-\\k-\  ((p),  P°(x),  z1:m] 

•  Predicted  target  intensity :  For  xeS, 

-  Newly  born  targets  (target  initiations):  bk(x)  =  ’fp-iMP)/jt-i|it-i(p) 

-  Moving  targets:  dk(x)  =  fs  %-i(xly)fk-ljk-1(i/)dy 

-  Target  intensity:  fk\k-i(x)  =  h(x)  +  dk{x) 

•  Predicted  clutter  intensity: 

-  Clutter  persistence:  bk(tp)  =  ’F|t-i(p|p)./jc-i|*:-i(p) 

-  Newly  born  clutter  (target  terminations):  dk(tp )  =  ¥4-i(p]y)/jt-i|it-i(y)dy 

-  Clutter  intensity:  fk\k-d<p)  =  h(<P)  +  dk(<p) 

•  Predicted  measurement  intensity:  FOR  j  =  l:m, 

-  Generated  by  targets:  vk{zj)  =  fs  pk{zj\x)P°(x)  fk[k_ i(x)dx 

-  Generated  by  clutter:  Xk(zj)  =  pk(Zj\<p)P°(cp)  fk\k-i((p) 

-  Measurement  intensity:  Ak\k-i(zj)  =  Ak{zj)  +  vk(zj) 

•  Information  updated  target  intensity :  For  xeS, 

-  Undetected  targets:  /^(x)  =  (l  -  Pk(x))fk\k-i(x) 

-  Detected  targets:  /°(x)  =  [l;m=1  lW 

-  Target  intensity:  fm(x)  =  /^(x)  +  /°  (x) 

•  Information  updated  clutter  intensity: 

-  Undetected  clutter  targets  (generate  no  data):  /^(p)  =  (l  -  Pj?(p))//c|it-i(p) 

-  Detected  clutter  targets  (generate  data):  f°k(<p)  =  A|fc-i (<P) 

-  Clutter  target  intensity:  _%(p)  =  /^(p)  +  /*>  (p) 


where 


Xk\k-l(z) 


PkizWpVWfw-iix)  dx, 


Ak{z)  =  pk(z\(p)Pk((p)fk\k-i((p) 


(6.15) 


(6.16) 


is  the  predicted  measurement  clutter  intensity.  The  probability  Pk(<p)  in  (6.16) 
is  the  probability  that  a  clutter  target  generates  a  measurement  at  time  tk- 
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The  computational  parts  of  the  intensity  filter  are  outlined  in  Table  6.1. 
The  table  clarifies  certain  interpretive  issues  that  are  glossed  over  in  the 
discussion.  Implementation  methods  are  discussed  elsewhere  in  this  chapter. 


Likelihood  Function  of  the  Data  Set 


Since  fm(x)  is  the  intensity  of  a  PPP,  it  is  reasonable  to  inquire  about  the 
pdf  of  the  data  r\k  =  (mk,  {zi,  • . .,  The  measurement  intensity  after  the 
information  update  of  the  target  state  intensity  is,  applying  (2.86)  on  the 
augmented  space  S+ , 


Pk(z\x)fk\k(x)dx. 


Therefore  the  pdf  of  the  data  r\k  is 


where 


mk 

P(r\k)  =  Ak*k{z)dz  Ak\k(zj) 


;= i 


mk 


=  e~N*  J]  Afc| k(zj), 

;=i 


Nk]k  =  I  Ajt|jt(z)dz 

=  ||  Pfc(z  |  x)  dz  ]  fk\k(x)  dx 

Js+  \Jr"z 

=  I  fm(x)dx 

J  Rnx 


(6.17) 


(6.18) 


(6.19) 


is  the  estimated  mean  number  of  targets.  The  pdf  of  rjjt  is  approximate  because 
fk\k(x)  is  an  approximation. 


6.2  Relationship  to  Other  Filters 

The  modeling  assumptions  of  the  intensity  filter  are  very  general,  and  spe¬ 
cializations  are  possible.  The  most  important  is  the  PHD  filter  discussed  in 
Section  6.2.1.  By  assuming  certain  kinds  of  a  priori  knowledge  concerning 
target  birth  and  measurement  clutter,  and  adjusting  the  filter  appropriately, 
the  intensity  filter  reduces  to  the  PHD  filter.  The  differences  between  the 
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intensity  and  PHD  filters  are  nearly  all  attributable  to  the  augmented  state 
space  »S+ .  That  is,  the  intensity  filter  uses  the  augmented  single  target  state 
space  *S+  =  S  U  cp,  while  the  PHD  filter  uses  only  the  single  target  space  S. 
Using  S  practically  forces  the  PHD  filter  to  employ  target  birth  and  death 
processes  to  model  initiation  and  termination  of  targets. 

A  different  kind  of  specialization  is  the  marked  multitarget  intensity  filter. 
This  is  a  parameterized  linear  Gaussian  sum  intensity  filter  that  interprets 
measurements  as  target  marks.  This  interpretation  is  interesting  in  the  con¬ 
text  of  PPP  target  models  because  it  implies  that  joint  measurement-target 
point  process  is  a  PPP.  Details  are  discussed  in  Section  6.2.2. 


6.2.1  Probability  Hypothesis  Density  (PHD)  Filter 

The  state  cp  is  the  basis  for  the  on-line  estimates  of  the  intensities  of  the  target 
birth  and  measurement  clutter  PPPs  given  by  (6.13)  and  6.15),  respectively. 
If,  however,  the  birth  and  clutter  intensities  are  known  a  priori  to  be  Iy(x)  and 
Ajt(z),  then  the  predictions  bk(x)  and  Ajt(z)  can  be  replaced  by  b] t(x)  and  Ajt(z). 
This  is  the  basic  strategy  taken  by  the  PHD  filter. 

The  use  of  a  posteriori  methods  makes  good  sense  in  many  applications. 
For  example,  they  can  help  regularize  parameter  estimates.  These  methods 
can  also  incorporate  information  not  included  a  priori  in  the  Bayes  filter.  For 
example,  Jazwinski  [49]  uses  an  a  posteriori  method  to  derive  the  Schmidt- 
Kalman  filter  for  bias  compensation.  These  methods  may  improve  perfor¬ 
mance,  i.e.,  if  the  a  priori  birth  and  clutter  intensities  are  more  accurate  or 
stable  than  their  on-line  estimated  counterparts,  the  PHD  filter  may  provide 
better  tracking  performance. 

Given  these  substitutions,  the  augmented  space  is  no  longer  needed  and 
can  be  eliminated.  This  requires  some  care.  If  the  recursion  is  simply  re¬ 
stricted  to  S  and  no  other  changes  are  made,  the  filter  will  not  be  able  to 
discard  targets  and  the  target  count  may  balloon  out  of  control.  To  balance 
the  target  birth  process,  the  PHD  filter  uses  a  death  probability  before  propa¬ 
gating  the  multitarget  intensity  fk-i\k-i(x).  This  probability  was  intentionally 
omitted  from  the  intensity  filter  because  transition  into  cp  is  target  death,  and 
it  is  redundant  to  have  two  death  models. 

The  death  process  is  a  Bernoulli  thinning  process  applied  to  the  PPP 
at  time  fj.- 1  before  targets  transition  and  are  possibly  detected.  Let  t4-i  (x) 
denote  the  probability  that  a  target  at  time  4_i  dies  before  transitioning 
to  time  tk .  The  surviving  target  point  process  is  a  PPP  and  its  intensity  is 
(1  -  dk-\(x))  fk-\\k-i(x).  Adding  Bernoulli  death  and  restricting  the  recursion 
to  S  reduces  the  intensity  filter  to  the  PHD  filter. 
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Table  6.2  PHD  Filter  on  the  State  Space  S 
INPUTS: 

Data:  z\:m  =  \z\,  c  T  at  time  tk 

Target  death  probability  function:  4-1  (x)  =  Pr  [target  death  at  state  x  e  S  at  time  f jt_i  ] 
Target  birth  probability  function:  bk(x)  =  Pr  [target  birth  at  state  x  €  S  at  time  tk] 
Probability  of  target  detection:  P°(x)  at  time  tk 
Measurement  clutter  intensity  function:  Ak(z)  at  time  tk 

OUTPUT:  fm{x)  =  PHDFilter (x),  4_ 1(*),  h(x),  Ak(z),  P°(x),  z1:m] 

•  Predicted  target  intensity :  For  xeS, 

-  Surviving  targets: 

Sk(x)  =  (i  -  4-iW)4-i|t-i(y) 

-  Propagated  targets: 


Sk(x)  «-  Js  Wk-i(x\y)Sk(y)dy 


-  Predicted  target  intensity: 

fk\k-i(x)  =  bk(x)  +  Sjt(x) 

•  Predicted  measurement  intensity: 

IF  m  =  0,  THEN  fm(x)  =  (l  -  Jf(x))  4,n(x)  STOP 
FOR  y'  =  l:  m, 

-  Intensity  contributions  from  predicted  target  intensity: 

n(Z;)  =  JsPk(zj\x)P^(x)fklk-1(x)dx 

-  Predicted  measurement  intensity: 

Ak\k-i(zj)  =  A  k(zj)  +  vk{zj) 

END  FOR 

•  Information  updated  target  intensity:  For  xeS, 


•  END 


6.2.2  Marked  Multisensor  Intensity  Filter  (MMIF) 

The  intensity  filter  assumes  targets  have  the  same  motion  model,  and  that 
the  sensor  measurement  likelihood  function  is  the  same  for  all  targets  and 
data.  Such  assumptions  are  idealized  at  best.  An  alternative  approach  is  to 
develop  a  parameterized  intensity  filter  that  accommodates  heterogeneous 
target  motion  models  and  measurement  pdfs  by  using  target-specific  param- 
eterizations.  The  notion  of  target-specific  parameterizations  in  the  context  of 
PPP  target  modeling  seems  inevitably  to  lead  to  the  idea  of  modeling  indi¬ 
vidual  targets  as  a  PPP,  and  then  using  superposition  to  obtain  the  aggregate 
PPP  target  model.  Parameter  estimation  using  the  EM  method  is  natural  to 
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superposition  problems,  as  shown  in  Chapter  3.  The  marked  multisensor 
intensity  filter  (MMIF)  is  one  instance  of  such  an  approach. 

The  MMIF  builds  on  the  basic  idea  that  a  target  at  state  x  is  "marked"  with 
a  measurement  z.  If  the  target  is  modeled  as  a  PPP,  then  the  joint  measurement- 
target  vector  (z,  x)  is  a  PPP  on  the  Cartesian  product  of  the  measurement  and 
target  spaces.  This  is  an  intuitively  reasonable  result,  but  the  details  needed 
to  see  that  it  is  true  are  postponed  to  Section  8.1.  The  MMIF  uses  a  lin¬ 
ear  Gaussian  target  motion  and  measurement  model  for  each  target  and 
superposes  them  against  a  background  clutter  model.  Since  the  Gaussian 
components  correspond  to  different  targets,  they  need  not  have  the  same 
motion  model.  Similarly,  different  sensor  measurement  models  are  possible. 
Superposition  therefore  leads  to  an  affine  Gaussian  sum  intensity  function 
on  the  joint  measurement-target  space.  The  details  of  the  EM  method  and 
the  final  MMIF  recursion  are  given  in  Appendix  E. 

The  MMIF  adheres  to  the  "at  most  one  measurement  per  target  rule"  but 
only  in  the  mean,  or  on  average.  It  does  this  by  reinterpreting  the  single 
target  pdf  as  a  PPP  intensity  function,  and  by  interpreting  measurements  as 
the  target  marks.  The  expected  number  of  targets  that  the  PPP  on  the  joint 
measurement-target  space  produces  is  one. 

Another  feature  of  the  MMIF  is  that  the  EM  weights  depend  on  the  Kalman 
filter  innovations.  The  weights  in  other  Gaussian  sum  filters  often  involve 
scaled  multiples  of  the  measurement  variances,  resulting  in  filters  that  are 
somewhat  akin  to  "nearest  neighbor"  tracking  filters. 

The  limitation  of  MMIF  and  other  parameterized  sum  approaches  is  the 
requirement  to  use  a  fixed  number  of  terms  in  the  sum.  This  strongly  affects 
its  ability  to  model  the  number  of  targets.  In  practice,  various  devices  can 
compensate  for  this  limitation,  but  they  are  not  intrinsic  to  the  filter. 


6.3  Implementation 

Simply  put,  targets  correspond  to  the  local  peaks  of  the  intensity  function 
and  the  areas  of  uncertainty  correspond  to  the  contours,  or  isopleths,  of  the 
intensity.  Very  often  in  practice,  isopleths  are  approximated  by  ellipsoids  in 
target  state  space  corresponding  to  error  covariance  matrices.  Methods  for 
locating  the  local  peak  concentrations  of  intensity  and  finding  appropriate 
covariance  matrices  to  measure  the  width  of  the  peaks  are  discussed  in  this 
section. 

Implementation  issues  for  intensity  filters  therefore  concern  two  issues. 
Firstly,  it  is  necessary  to  develop  a  computationally  viable  representation  of 
the  information  updated  intensity  function  of  the  filter.  Two  basic  representa¬ 
tions  are  proposed,  one  based  on  particles  and  the  other  on  Gaussian  sums. 
Secondly,  postprocessing  procedures  are  applied  to  the  intensity  function 
representation  to  extract  the  number  of  detected  targets,  together  with  their 
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estimated  states  and  corresponding  error  covariance  matrices.  Analogous 
versions  of  both  issues  arise  in  classical  single  target  Bayesian  filters. 

The  fact  remains,  however,  that  a  proper  statistical  interpretation  of  target 
point  estimates  and  their  putative  error  covariances  is  lacking  for  intensity 
filters.  The  concern  may  be  dismissed  in  practice  because  they  are  intuitively 
meaningful  and  closely  resemble  their  single  target  Bayesian  analogs.  The 
concern  is  nonetheless  worrisome  and  merits  further  study. 


6.3.1  Particle  Methods 

The  most  common  and  by  far  the  easiest  implementation  of  nonlinear  filters 
is  by  particle,  or  sequential  Monte  Carlo  (SMC),  methods.  In  such  methods 
the  posterior  pdf  is  represented  nonparametrically  by  a  set  of  particles  in 
target  state  space,  together  with  a  set  of  associated  weights,  and  estimated 
target  count.  Typically  these  weights  are  uniform,  so  the  spatial  distribution 
of  particles  represents  the  variability  of  the  posterior  density.  An  excellent 
discussion  of  SMC  methods  for  Bayesian  single  target  tracking  applications 
is  found  in  the  first  four  chapters  of  [91]. 

Published  particle  methods  for  the  general  intensity  filter  are  limited  to 
date  to  the  PHD  filter.  Extensions  to  the  intensity  filter  are  not  reported  here. 
An  early  and  well  described  particle  methodology  (as  well  as  an  interesting 
example  for  tracking  on  roads)  for  PHD  filters  is  given  in  [98].  Particle 
methods  and  their  convergence  properties  for  the  PHD  filter  are  discussed 
in  detail  in  a  series  of  papers  by  Vo  et  al.  [126].  Interested  readers  are  urged 
to  consult  them  for  specifics. 

Tracking  in  a  surveillance  region  R  using  SMC  methods  starts  with  an 
initial  set  of  particles  and  weights  at  time  4_i  together  with  the  estimated 
number  of  targets  in  R: 


where  \\k~i(t)  =  1  /Rsmc  for  all  L  For  PHD  filters  the  particle  method 
proceeds  in  several  steps  that  mimic  the  procedure  outlined  in  Table  6.2: 

•  Prediction.  In  the  sequential  importance  resampling  (SIR)  method,  predic¬ 
tion  involves  thinning  a  given  set  of  particles  with  the  survival  probability 
1  -  dRx)  and  then  stochastically  transforming  the  survivors  into  another 
particle  set  using  the  target  motion  model  V'jt- i(j/|x)  with  weights  ad¬ 
justed  accordingly.  Additional  new  particles  and  associated  weights  are 
generated  to  model  new  target  initializations. 

•  Updating.  The  particle  weights  are  multiplicatively  (Bayesian)  updated 
using  the  measurement  likelihood  function  p^(z  \  x)  and  the  probability  of 
detection  P^(x).  The  factors  are  of  the  form 


5-18 


RTO-EN-SET-1 57(2010) 


Poisson  Point  Processes:  Imaging,  Tracking,  and  Sensing 


168 


6  Applications:  Multiple  Target  Tracking 


(6.20) 


where  zjt(l), . . . ,  z /,(»/)  are  the  measurements  at  time  />.  The  updated  par¬ 
ticle  weights  are  nonuniform. 

•  Normalization.  Compute  the  scale  factor,  call  it  Nk\k,  of  the  sum  of  the  up¬ 
dated  particle  weights.  Divide  all  the  particle  weights  by  Nk\k  to  normalize 
the  weights. 

•  Resampling.  Particles  are  resampled  by  choosing  i.i.d.  samples  from  the 
discrete  pdf  defined  by  the  normalized  weights.  Resampling  restores  the 
particle  weights  to  uniformity. 

If  the  resampling  step  is  omitted,  the  SMC  method  leads  to  particle  weight 
distributions  that  rapidly  concentrate  on  a  small  handful  of  particles  and 
therefore  poorly  represents  the  posterior  intensity.  There  are  many  ways  to 
do  the  resampling  in  practice. 

By  computing  Nk\k  before  resampling,  it  is  easy  to  see  that 


Jk 

=  E  [Number  of  targets  in  R\  . 

The  estimated  number  of  targets  in  any  given  subset  Rq  C  R  is 


(6.21) 


(6.22) 


The  estimator  is  poor  for  sets  Ro  that  are  only  a  small  fraction  of  the  total 
volume  of  R. 

The  primary  limitations  of  particle  approaches  in  many  applications  are 
due  to  the  so-called  Curse  of  Dimensionality1 :  the  number  of  particles  needed 
to  represent  the  intensity  function  grows  exponentially  as  the  dimension 
of  the  state  space  increases.  Most  applications  to  date  seem  to  be  limited 
to  four  or  five  dimensions.  The  curse  is  so  wicked  that  Moore's  Law  (the 
doubling  of  computational  capability  every  18  months)  by  itself  will  do  little 
to  increase  the  effective  dimensional  limit  over  a  human  lifetime.  Moore's 
Law  and  improved  methods  together  will  undoubtedly  increase  the  number 
of  dimensions  for  which  particle  filters  are  practical,  but  it  remains  to  be 
seen  if  general  filters  of  dimension  much  larger  than  say  six  can  be  treated 
directly. 


1  The  name  was  first  used  in  1961  by  Richard  E.  Bellman  [8].  The  name  is  apt  in  very 
many  problems;  however,  some  modern  methods  in  machine  learning  actually  exploit 
high  dimensional  embeddings. 
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6.3.2  Mean  Shift  Algorithm 

As  is  clear  from  the  earlier  discussion  of  PET,  it  is  intuitively  reasonable 
to  think  of  the  multitarget  intensity  filter  as  a  sequential  image  processing 
method.  In  this  interpretation  an  image  comprises  the  “gray  scales"  of  a  set  of 
multidimensional  voxels  in  the  target  state  space  R”x .  Such  an  interpretation 
enables  a  host  of  image  processing  techniques  to  be  applied  to  the  estimated 
intensity  function. 

One  technique  that  lends  itself  immediately  to  extracting  target  point 
estimates  from  a  particle  set  approach  is  the  "mean-shift"  algorithm.  This 
algorithm,  based  on  ideas  first  proposed  in  [35],  is  widely  used  in  computer 
vision  applications  such  as  image  segmentation  and  tracking.  The  mean-shift 
algorithm  is  an  EM  algorithm  for  Gaussian  kernel  density  estimators  and 
a  generalized  EM  algorithm  for  non-Gaussian  kernels  [11].  The  Gaussian 
kernel  is  computationally  very  efficient  in  the  mean-shift  method. 

Denote  the  non-([)  particle  set  representing  the  PPP  intensity  at  time  4  by 


The  intensity  function  is  modeled  as  a  scalar  multiple  of  the  kernel  estimator 


(6.23) 


where  ;  Xk\k{t),  Hker)  is  the  kernel.  The  covariance  matrix  I’ker  is  specified, 

not  estimated.  Intuitively,  the  larger  Zker/  the  fewer  the  number  of  local 
maxima  in  the  intensity  (6.23),  and  conversely.  The  scale  factor  4  >  0  is 
estimated  by  the  particle  filter  and  is  taken  as  known  here. 

The  form  (6.23)  has  no  parameters  to  estimate,  so  extend  it  by  defining 


(6.24) 


where  p is  an  unknown  rigid  translation  of  the  intensity  (6.23).  It  is  not  hard 
to  see  that  the  ML  estimate  of  p  is  a  local  maximum  of  the  kernel  estimate, 
that  is,  a  point  estimate  for  a  target.  The  vector  p  is  estimated  from  data  using 
the  EM  method.  The  clever  part  is  using  an  artificial  data  set  with  only  one 
point  in  it,  namely,  the  origin. 

Let  r  =  0, 1, . . .  denote  the  EM  iteration  index,  and  let  p((l)  be  a  specified 
initial  value  for  the  mean.  The  auxiliary  function  is  given  by  Eqn.  (3.20)  with 
m  =  1,  xx  =  0,  L  =  Lsmc,  dc  =  p,  and 
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The  bounded  surveillance  region  R  is  taken  to  be  1R"1.  Define  the  weights 

^(0;xklM-^\LkeT) 

iff  ;  =  — - - - - - 

v  ’  L^c^(0;xk{k(n-^,^Cr) 

=  - 7 - 7  •  (6-25) 

The  auxiliary  function  in  the  present  case  requires  no  sum  over  j  as  done  in 
(3.20),  so 


lsmc 


Q(fb  p(r))  =  -Ik+  Wf(/-‘(r))  log {4^(0;xfc|jt(^)  -  ft,  £ker)} .  (6.26) 


t= l 


The  EM  update  of  ft  is  found  by  taking  the  appropriate  gradient,  yielding 

^(fi(r))  xwd) 


f' 


(r+l)  _ 


L^C  ^(pW) 


(6.27) 


Substituting  (6.25)  and  canceling  the  common  factor  gives  the  classical  mean- 
shift  iteration: 


fr+1 ,  LlAT  ^  (xkVcW ;  ft(r),  2iker)  xm{{) 

ft1  =  - - 7 - 7 - •  (6.28) 

LLlT^(xMM-r^\Ekei) 

The  update  of  the  mean  is  a  convex  combination  of  the  particle  set.  Conver- 

(r) 

gence  to  a  local  maximum  pjj.  — >  xkk  is  guaranteed  as  r  — >  oo. 

Different  initializations  are  needed  for  different  targets,  so  the  mean  shift 
algorithm  needs  a  preliminary  clustering  method  to  initialize  it,  as  well  as  to 
determine  the  number  of  peaks  in  the  data  correspond  to  targets.  Also,  the 
size  of  the  kernel  depends  somewhat  on  the  number  of  particles  and  may 
need  to  be  adjusted  to  smooth  the  intensity  surface  appropriately. 


6.3.3  Multimode  Algorithms 

Identifiability  remains  a  problem  with  the  mean  shift  algorithm,  that  is,  there 
is  no  identification  of  the  point  estimate  to  a  target  except  through  the  way 
the  starting  point  of  the  iteration  is  chosen.  This  may  cause  problems  when 
targets  are  in  close  proximity. 

One  way  to  try  to  resolve  the  problem  is  to  use  the  particles  themselves 
as  points  to  feed  into  another  tracking  algorithm.  This  method  exploits  se- 
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rial  structure  in  the  filter  estimates  and  may  disambiguate  closely  spaced 
targets.  However,  particles  are  serially  correlated  and  do  not  satisfy  the  con¬ 
ditional  independence  assumptions  of  measurement  data,  so  the  resulting 
track  estimates  may  be  biased. 


6.3.4  Covariance  Matrices 


Matrix  CRBs 


The  true  error  covariance  matrix  of  the  ML  point  estimate  computed  by  the 
mean  shift  method  is  not  available;  however,  the  CRB  can  be  evaluated  using 
(4.21).  This  CRB  is  appropriate  if  the  available  data  are  reasonably  modeled 
as  realizations  of  a  PPP  with  intensity  (6.24).  If  the  intensity  (6.24)  function 
is  a  high  fidelity  model  of  the  intensity,  the  FIM  of  /./  is 


m  =  h 


(6.29) 


where  the  matrix  E(f.i)  is 


L  L  r  T 

E(j.t)  =  V  V  I  zv(x>  ;  C  C)[x- xk\k(t)  +  fij (x - xk\k({')  +  j.i'j  dx  (6.30) 

f=i  c= l  ^]Rnx 


and  the  weighting  function  is 


iv (x,  p  ■,(,{') 


{X',  Xk\k(i)  ^ker) 

\jt\ 

[x;xm({')  - 

/ -b  -^ker) 

^X)X 

k\k(t")  ~  b/  A;er) 

i 

(6.31) 


The  CRB  is  J_1(p)  evaluated  at  the  true  value  of  p.  Because  the  integral  is 
over  all  of  IR"1,  a  change  of  variables  shows  that  the  information  matrix  /(p) 
is  independent  of  the  true  value  of  p.  This  means  that  the  FIM  is  not  target 
specific. 

A  local  bound  is  desired,  since  the  mean  shift  algorithm  converges  to  a 
local  peak  of  the  intensity.  By  restricting  the  intensity  model  to  a  specified 
bounded  gate  G  C  IR"’:,  the  integral  in  (6.30)  is  similarly  restricted.  The  matrix 
£(p)  is  thus  a  function  of  G.  The  gated  CRB  is  local  to  the  gate,  i.e.,  it  is  a 
function  of  the  target  within  the  gate. 


OIM:  The  Surrogate  CRB 

The  OIM  is  the  Hessian  matrix  of  the  negative  loglikelihood  evaluated  at  the 
MAP  point  estimate;  it  is  often  used  as  a  surrogate  for  Fisher  information 
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when  the  likelihood  function  is  complicated.  Its  inverse  is  the  surrogate  CRB 
of  the  estimate  x^.  The  construction  of  OIMs  for  the  mean-shift  algorithm  is 
implicit  in  [11]. 

The  loglikelihood  function  of  fi  using  the  intensity  function  (6.23)  and  the 
one  data  point  X\  =  0  is 


log  p(fO  =  -IkLsMC  +  log 


t-SMC 

^  ^  ‘  1  (Xk\k(0  /  y /  Tker) 
e=i 


Direct  calculation  gives  the  general  expression 

i  \Lsmc 


Vfi[Vfilog  P(fl)]  =£; 


ker  +  k2 


Xk\k(C),  Tker)TkfIr  (ft  -  Xty(ej) 


t= 1 


lsmc 


_ i;-1 


ker 


^ (ft ;  xm({),  Zker)  hjr  (ft  -  x*|fc(£)) 
f=i 
lsmc 

^  ^  '  1  (fl  /  -Tk|k(f)/  Tker)  (ft  —  Xk|k(t)  j  (f/  —  Vk|k(t)  j 


f=l 


y-l 

^ker' 


where  the  normalizing  constant  is 

lsmc 


K  = 


■  T(ft  /  xk|k(f),  Tker) . 


£=1 


The  observed  information  matrix  is  evaluated  at  the  MAP  estimate  ft  =  Xm. 

The  middle  term  is  proportional  to  (Vf,  p(ft)j  (Vf,  p(ft)j  and  so  is  zero  any 
stationary  point  of  p(ft),  e.g.,  at  the  MAP  estimate  ft  =  xk|k.  The  OIM  is 
therefore 


OIM(xm)  = 

%k\k(,^)/  ^ker)  (%k\k  ~  Xk\ki^))  (%k\k  — 


-  r_1 

^ker 


LtlT  ker) 


y— 1 

^“'ker  * 


(6.32) 


The  CRB  surrogate  is  OfM  1  (xk|k).  The  inverse  exists  because  the  OIM  is 
positive  definite  at  Xk\k-  This  matrix  is,  in  turn,  a  surrogate  for  the  error 
covariance  matrix. 

The  OIM  for  Xk\k  can  be  computed  efficiently  in  conjunction  with  any 
EM  method  (see  [61]  for  a  general  discussion).  As  noted  in  Section  4.7,  the 
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statistical  interpretation  of  the  OIM  is  unresolved  in  statistical  circles.  Its 
utility  should  be  carefully  investigated  in  applications. 


6.3.5  Gaussian  Sum  Methods 

An  alternative  to  representing  the  PPP  intensity  using  particle  filters  is  to 
use  a  Gaussian  sum  instead.  An  advantage  of  this  method  is  that  the  clus¬ 
tering  and  track  extraction  issues  are  somewhat  simpler,  that  is,  target  state 
estimates  and  covariance  matrices  are  extracted  from  the  means  and  vari¬ 
ances  of  the  Gaussian  sum  instead  of  a  myriad  of  particles.  The  Gaussian 
sum  approach  is  especially  attractive  for  linear  Gaussian  target  motion  and 
measurement  models  because  the  prediction  and  information  update  steps 
are  closed  form,  assuming  constant  survival  and  detection  functions  (i.e.,  the 
PPP  thinning  functions  are  independent  of  state). 

Gaussian  sum  implementations  of  the  PHD  filter  are  carefully  discussed 
by  Vo  and  his  colleagues  [124].  In  this  approach,  an  unnormalized  Gaussian 
sum  is  used  to  approximate  the  intensity.  These  methods  are  important 
because  they  have  the  potential  to  be  useful  in  higher  dimensions  than 
particle  filters.  In  the  end,  however,  Gaussian  sum  intensity  filters  will  also 
suffer  from  the  curse  of  target  state  space  dimensionality. 

Gaussian  sum  methods  for  intensity  estimation  comprises  several  steps: 

•  Prediction.  The  target  intensity  at  time  fj.- 1  is  a  Gaussian  sum,  to  which 
is  added  a  target  birth  process  that  is  modeled  by  a  Gaussian  sum.  The 
prediction  equation  for  every  component  in  the  Gaussian  sum  is  identical 
to  a  Kalman  filter  prediction  equation. 

•  Component  Update.  For  each  point  measurement,  the  predicted  Gaussian 
components  are  updated  using  the  usual  Kalman  update  equations.  The 
update  therefore  increases  the  number  of  terms  in  the  Gaussian  sum  if 
there  is  more  than  one  measurement.  This  step  has  two  parts.  In  the  first, 
the  means  and  covariance  matrices  are  evaluated.  In  the  second,  the  coef¬ 
ficients  of  the  Gaussian  sum  are  updated  by  a  multiplicative  procedure. 

•  Merging  and  Pruning.  The  components  of  the  Gaussian  sum  are  merged 
and  pruned  to  obtain  a  "nominal"  number  of  terms.  Various  reasonable 
strategies  are  available  for  such  purposes,  as  detailed  in  [124].  This  step  is 
the  analog  of  resampling  in  the  particle  method. 

Some  form  of  pruning  is  necessary  to  keep  the  size  of  the  Gaussian  sum 
bounded  over  time,  so  the  last  -  and  most  heuristic  -  step  cannot  be  omitted. 

Left  out  of  this  discussion  are  details  that  relate  the  weights  of  the  Gaussian 
components  to  the  estimated  target  count.  These  details  can  be  found  in  [124]. 
For  nonlinear  target  motion  and  measurement  models,  [124]  proposes  both 
the  extended  and  the  unscented  Kalman  filters.  Vo  and  his  colleagues  also 
present  Gaussian  sum  implementations  of  the  CPHD  filter  in  [124]  and  [125]. 
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6.3.6  Regularization 

Intensity  filters  are  in  the  same  class  of  stochastic  inverse  problems  as  image 
reconstruction  in  emission  tomography  —  the  sequence  to,  h, . . . ,  tk  of  inten¬ 
sity  filter  estimates  fk\k(x)  is  essentially  a  movie  (in  dimension  nx)  in  target 
state  space.  As  discussed  in  Section  5.7,  such  problems  suffer  from  serious 
noise  and  numerical  artifacts.  The  high  dimensionality  of  the  PPP  parameter, 
i.e.,  the  number  of  voxels  of  the  intensity  function,  makes  regularization  a 
priority  in  all  applications.  Regularization  for  intensity  filters  is  a  relatively 
new  subject.  Methods  such  as  cardinalization  are  inherently  regularizing. 

Grenander's  method  of  sieves  used  in  Section  5.7  for  regularizing  PET 
adapts  to  the  intensity  filter,  but  requires  some  additional  structure.  The 
sieve  kernel  ko(x\u)  is  a  pdf  on  <S+,  so  that 


I  ko(x\u)dx  =  1 

Js+ 


(6.33) 


for  all  points  u  in  the  discrete-continuous  space  rU'  =  'll  U  (p.  As  before,  the 
choice  of  kernel  and  space  T/+  is  very  flexible.  The  multitarget  intensity  at 
every  time  tk,  k  =  0, 1, ...,  is  restricted  to  the  collection  of  functions  of  the 
form 


fk\k(x)  =  ko(x\u)  ik{k(u)du  for  some  L,k\k(u)  >  0.  (6.34) 

JU+ 


The  kernel  /cq  can  be  a  function  of  time  f;c  if  desired.  The  restriction  (6.34)  is 
also  imposed  on  the  predicted  target  intensity: 


fk\k-\(x)  ~  ...  ko(x|w)Cjt[Jc-i(w)dit  for  some  Cit|i:-i(u)  >  0-  (6.35) 


Substituting  (6.35)  into  the  predicted  measurement  intensity  (6.3)  gives 


(6.36) 


where 


(6.37) 


is  the  regularized  measurement  likelihood  function. 

An  intensity  filter  is  used  to  update  Citi k(u).  This  filter  employs  a  transition 
function  <T>fc_i(  ■  |  ■ )  to  provide  a  dynamic  connection  between  the  current 
intensity  Cfc-i|fc-i(w)  and  predicted  intensity  L,k\k-i(u)-  The  function  CE>jt_i(  - 1  ■ ) 
is  specified  for  all  u  and  v  in  rU+  and  is,  in  principle,  any  reasonable  function, 
but  in  practice  is  linked  to  the  target  motion  model  T/jL_i( '  I '  )■ 
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One  way  to  define  <J>j.-i('  I  •)  requires  defining  an  additional  kernel.  Sub¬ 
stitute  (6.34)  into  the  predicted  detected  target  intensity  to  obtain 

fk\k-i(x)  =  I  fk-i\k-i(yWk-i(x\y)dy 

JS+ 

= Hi  koiy  I U)  Cfc— i|fc— 1  (u)  du  j  Wk^(x |  y)  d y 


=  f 


^/c-i^hOCfc-iDc-i^Od  U, 


where 


Wk-i(x\u)  =  f  xPk-i(xly)k0(ylu)dy. 
JS+ 

Define  a  Bayesian  kernel  ki(v\x)  so  that 

ki(v\x)dv  =  1 


f 


(6.38) 


(6.39) 


for  all  points  x  e  S+.  Like  the  sieve  kernel  ko(-),  the  Bayesian  kernel  /C|  (z>  |  x) 
is  very  flexible.  It  is  easily  verified  that  the  function 


Ojt-ifalu)  =  f  k1(v\x)Wk_1(x\u)dx 
JS+ 


(6.40) 


is  a  valid  transition  function  for  all  kg(  ■ )  and  k\  ( ■ ). 

Given  the  intensity  C*— lilt— i(m)  from  time  the  predicted  intensity 
t)t|fc-i00  at  time  tk  is  defined  by 

Cfc|fc-i(w)  =  ®fc-l(M|»)Qt-l|fc-l(P)d0. 

The  information  updated  intensity  Cit|it(u)  is  evaluated  via  the  intensity  filter 
using  the  regularized  measurement  pdf  (6.37)  and  the  predicted  measure¬ 
ment  intensity  (6.36).  The  regularized  target  state  intensity  at  time  tk  is  the 
integral 


fk\k(x)  =  k(x\n)lklk(u)du.  (6.41) 

J-W+ 

The  regularized  intensity  fk\k(x)  depends  on  the  sieve  and  Bayesian  kernels 
fco(-)  and  MO- 

The  question  of  how  best  to  define  the  /cq(  • )  and  k\ ( ■ )  kernels  depends 
on  the  application.  It  is  common  practice  to  define  kernels  using  Gaussian 
pdfs.  As  mentioned  in  Section  5.7,  the  sieve  kernel  is  a  kind  of  measurement 
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smoothing  kernel.  If  dim('W)  <  dim(<S),  the  Bayesian  kernel  disguises  ob¬ 
servability  issues,  that  is,  many  points  x  e  S  map  with  the  same  probability 
to  a  given  point  u  e  14.  This  provides  a  mechanism  for  target  state  space 
smoothing. 


6.4  Estimated  Target  Count 

Estimating  the  number  of  targets  present  in  the  data  is  a  difficult  hypothesis 
testing  problem.  In  a  track  before  detect  (TBD)  approach,  it  is  a  track  man¬ 
agement  function  that  is  integrated  with  the  tracking  algorithm.  Intensity 
filters  seem  to  offer  an  alternative  way  to  integrate,  or  fuse,  the  multitarget 
track  management  and  state  estimation  functions. 


6.4.2  Sources  of  Error 

Accurate  knowledge  of  the  target  detection  probability  function  P^(x)  is 
crucial  to  correctly  estimating  target  count.  An  incorrect  value  of  Pjr(x)  is  a 
source  of  systematic  error.  For  example,  if  the  filter  uses  the  value  Pf(x)  =  .5 
but  in  fact  nil  targets  alzvays  show  up  in  the  measured  data,  the  estimated 
mean  target  count  will  be  high  by  a  factor  of  two.  This  example  is  somewhat 
extreme  but  it  makes  the  point  that  correctly  setting  the  detection  probability 
is  an  important  task  for  the  track  management  function.  The  task  involves 
executive  knowledge  about  changing  sensor  performance  characteristics, 
as  well  as  executive  decisions  external  to  the  tracking  algorithm  about  the 
number  of  targets  actually  present  —  decisions  that  feedback  to  validate 
estimates  of  P®(x).  Henceforth,  the  probability  of  detection  function  Pj?(x)  is 
assumed  accurate. 

There  are  other  possible  sources  of  error  in  target  count  estimates.  Birth- 
death  processes  can  be  difficult  to  tune  in  practice,  regardless  of  whether  they 
are  modeled  implicitly  as  transitions  into  and  out  of  a  state  <ft  in  the  intensity 
filter,  or  explicitly  as  in  the  PHD  filter.  If  in  an  effort  to  detect  new  targets 
early  and  hold  track  on  them  as  long  as  possible,  births  are  too  spontaneous 
and  deaths  are  too  infrequent,  the  target  count  will  be  too  high  on  average. 
Conversely,  it  will  be  too  low  with  delayed  initiation  and  early  termination. 
Critically  damped  designs,  however  that  concept  is  properly  defined  in  this 
context,  would  seem  desirable  in  practice.  In  any  event,  tuning  is  a  function 
of  the  track  management  system. 

Under  the  PPP  model,  the  estimated  expected  number  of  targets  in  a 
given  region,  J{,  is  the  integral  over  H  of  the  estimated  multitarget  intensity. 
Because  the  number  is  Poisson  distributed,  the  variance  of  the  estimated 
number  of  targets  is  equal  to  the  mean  number.  This  large  variance  is  an 
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unhappy  fact  of  life.  For  example,  if  10  targets  are  present,  the  standard 
deviation  on  the  estimated  number  is  VlO  ~  3 . 

It  is  therefore  foolhardy  in  practice  to  assume  that  the  estimated  of  the 
number  of  targets  is  the  number  that  are  actually  present.  Variance  reduction 
in  the  target  count  estimate  is  a  high  priority  from  the  track  management 
point  of  view  for  both  intensity  and  PHD  filters. 


6.4.2  Variance  Reduction 

The  multisensor  intensity  filter  discussed  in  Section  6.5  reduces  the  variance 
by  averaging  the  sensor-level  intensity  functions  over  the  number  of  sensors 
that  contribute  to  the  filter.  Consequently,  if  the  individual  sensors  estimate 
target  count  correctly,  so  does  the  multisensor  intensity  filter. 

Moreover,  and  just  as  importantly,  the  variance  of  the  target  count  estimate 
of  the  multisensor  intensity  filter  is  reduced  by  a  factor  of  M  compared  to  that 
of  a  single  sensor,  where  M  is  the  number  of  sensors,  assuming  for  simplicity 
that  the  sensor  variances  are  identical. 

This  important  variance  reduction  property  is  analogous  to  estimators 
in  other  applications.  An  especially  prominent  and  well  known  example 
is  power  spectral  estimation  of  wideband  stationary  time  series.  For  such 
signals  the  output  bins  of  the  DFT  of  a  non-overlapped  blocks  of  sampled 
data  are  distributed  with  a  mean  level  equal  to  the  signal  power  in  the 
bin,  and  the  variance  equal  to  the  mean.  This  property  of  the  periodogram 
is  well  known,  as  is  the  idea  of  time  averaging  the  periodogram,  i.e.,  the 
non-overlapped  DFT  outputs,  to  reduce  the  variance  of  spectral  estimates.2 
The  Wiener-Khinchin  theorem  justifies  averaging  the  short  term  Fourier 
transforms  of  nonoverlapped  data  records  as  a  way  to  estimate  the  power 
spectrum  with  reduced  variance.  In  practice,  the  number  of  DFT  records 
averaged  is  often  about  25. 

The  multisensor  intensity  filter  is  low  computational  complexity,  and  ap¬ 
plicable  to  distributed  heterogeneous  sensor  networks.  It  is  thus  practical  and 
widely  useful.  Speculating  now  for  the  sheer  fun  of  it,  if  the  number  of  data 
records  in  a  power  spectral  average  carries  over  to  multisensor  multitarget 
tracking  problems,  then  the  multisensor  intensity  filter  achieves  satisfactory 
performance  for  many  practical  purposes  with  about  25  sensors. 


2  Averaging  trades  off  variance  reduction  and  spectral  resolution.  It  was  first  proposed  by 
M.  S.  Bartlett  [6]  in  1948. 
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6.5  Multiple  Sensor  Intensity  Filters 

To  motivate  the  discussion,  consider  the  SPECT  application  of  Section  5.4.  In 
SPECT,  a  single  gamma  camera  is  moved  to  several  view  angles  and  a  snap¬ 
shot  is  taken  of  light  observed  emanating  from  gamma  photon  absorption 
events.  The  EM  recursion  given  by  Eqn.  (5.67)  is  the  superposition  of  the 
intensity  functions  estimated  by  each  of  the  camera  view  angles.  Intuitively, 
since  different  snapshots  cannot  contain  data  from  the  same  absorption  event, 
the  natural  way  to  fuse  the  multiple  images  into  one  image  is  to  add  them. 
The  theoretical  justification  is  that,  since  the  number  of  absorptions  is  un¬ 
known  and  Poisson  distributed,  the  estimates  of  the  spatial  distribution  of 
the  radioisotope  that  are  obtained  from  different  view  angles  are  indepen¬ 
dent,  not  conditionally  independent,  so  the  intensity  functions  (images)  are 
superposed. 

The  general  multisensor  filtering  problem  is  not  concerned  with  gamma 
photon  absorptions,  but  rather  with  physical  entities  (aircraft,  ships,  etc.)  that 
persist  over  long  periods  of  time  —  target  physics  are  very  different  from 
the  physics  of  photons.  Nonetheless,  for  reasons  discussed  at  the  beginning 
of  this  chapter,  a  PPP  multitarget  model  is  used  for  single  sensor  multitarget 
tracking.  It  is  assumed,  very  reasonably,  that  the  same  PPP  target  model  holds 
regardless  of  how  many  sensors  are  employed  to  detect  and  track  targets. 
The  analogy  with  SPECT  is  now  clear:  each  sensor  in  the  multisensor  filtering 
problem  is  analogous  to  a  camera  view  angle  in  SPECT,  and  the  sensor-level 
data  are  analogous  to  the  camera  snapshot  data. 

The  PPP  multitarget  model  has  immediate  consequences.  The  most  im¬ 
portant  is  that  conditionally  independent  sensors  are  actually  independent 
due  to  the  independence  property  of  PPPs  discussed  in  Section  2.9.  The 
multisensor  intensity  filter  averages  the  sensor-level  intensities  [111].  The 
reason  the  sensor  intensities  are  averaged  and  not  simply  added  is  that  tar¬ 
gets  are  persistent  (unlike  absorption  events).  Adding  the  intensities  "over 
count"  targets  because  each  sensor  provides  its  own  independent  estimate 
of  target  intensity. 

The  possibility  of  regularization  at  the  multisensor  level  is  not  considered 
explicitly.  Although  perhaps  obvious,  the  multisensor  intensity  filter  is  fully 
compatible  with  sensor-level  regularization  methods. 

Let  the  number  of  sensors  be  M  >  1 .  It  is  assumed  that  the  target  detection 
probability  functions,  P^(x;  £),  £  =  1, . . . ,  M,  are  specified  for  each  sensor.  The 
sensor-specific  state  space  coverage  is  defined  by 

Ck(£)  =  {x  e  S  :  P°(x ;  £)  >  o} .  (6.42) 

In  homogeneous  problems  the  sensor  coverages  are  identical,  i.e.,  Ck(£)  =  Ck 
for  all  L  Heterogeneous  problems  are  those  that  are  not  homogeneous. 

Two  sensors  with  the  same  coverage  need  not  have  the  same,  or  even 
closely  related,  probability  of  detection  functions.  As  time  passes,  homoge- 


RTO-EN-SET-1 57(2010) 


5-29 


Poisson  Point  Processes:  Imaging,  Tracking,  and  Sensing 


ORGANIZATION 


6.5  Multiple  Sensor  Intensity  Filters  179 

neous  problems  may  turn  into  heterogeneous  ones,  and  vice  versa.  In  prac¬ 
tice,  it  is  probably  desirable  to  set  a  small  threshold  to  avoid  issues  with  very 
small  probabilities  of  detection.  Homogeneous  and  heterogeneous  problems 
are  discussed  separately. 


6.5.1  Identical  Coverage  Sensors 


For  £  =  1, . . . ,  M,  let  the  measurement  space  of  sensor  (  be  £(■£).  Denote  the 
measurement  pdf  by  pk(z\x;{),  where  z{()  e  'Z.il)  is  a  point  measurement. 
The  predicted  measurement  intensity  is 

hk\k-i(z;  t)=  f  Pk(z\x;t)P°(x;{)fk |fc_i(x)dx,  z  e  Z(f)  ■  (6.43) 

Js+ 

The  measured  data  from  sensor  t  is 


4(f)  =  {zk(l  zk{mk(t)  ;  f )}) . 

The  sensor-level  intensity  filter  is 

fk\k(x>C)  =  4(4(f)l*;  O/fcifc-iW, 

where  the  Bayesian  information  update  factor  is 

n  pk{zk(i ;  €)  | x ;  t) Pf(x ;  €) 

The  multisensor  intensity  filter  is  the  average: 

M 


f=i 

M 


=  ^  E  4(4(f)U;f) 

u=i 


fk\k-l(x)- 


(6.44) 


(6.45) 


(6.46) 


(6.47) 


If  the  sensor-level  intensity  filters  are  maintained  by  particles,  and  the  num¬ 
ber  of  particles  is  the  same  for  all  sensors,  the  multisensor  averaging  filter  is 
implemented  merely  by  pooling  all  the  particles  (and  randomly  downsam¬ 
pling  to  the  desired  particle  size,  if  desired). 

Multisensor  fusion  methods  sometimes  rank  sensors  by  some  relative 
quality  measure.  This  is  unnecessary  for  the  multisensor  intensity  filter.  The 
reason  is  that  sensor  quality,  as  measured  by  the  probability  of  detection 
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functions  P;D(x;  (')  and  the  sensor  measurement  pdfs  pk(z  \  x ;  £),  is  automati¬ 
cally  included  in  (6.47). 

The  multisensor  intensity  filter  estimates  the  number  of  targets  as 

K?ed=  teseV)d* 

us 

i  M  r 

=  mH  fw(xU)dx 

1  e=i  ds 

=  (6-48) 

t=\ 

where  Nm{£)  is  the  number  of  targets  estimated  by  sensor  t.  Taking  the 
expectation  of  both  sides  gives 

1  M 

£Kr1]  =  ME£Kw]-  <6-49> 

t=\ 


If  the  individual  sensors  are  unbiased  on  average,  or  in  the  mean,  then 
E[Nk\k(()]  =  N  for  all  f,  where  N  is  the  true  number  of  targets  present.  Con¬ 
sequently,  the  multisensor  intensity  filter  is  also  unbiased. 

The  estimate  Nk\k(£)  is  Poisson  distributed,  and  the  variance  of  a  Poisson 
distribution  is  equal  to  its  mean,  so 


Var  [Nm({)]=N,  €  =  1  ,...,M. 


The  variance  of  the  average  in  (6.48)  is  the  average  of  the  variances,  since  the 
terms  in  the  sum  are  independent.  Thus, 


1  M 

VarKr4]  =  „  E 

f=l 

-  E 

_  M  ’ 


(6.50) 


In  words,  the  standard  deviation  of  the  estimated  target  count  in  the  multi¬ 
sensor  intensity  filter  is  smaller  than  that  of  individual  sensors  by  a  factor  of 
VM,  where  M  is  the  number  of  fielded  sensors.  This  is  an  important  result 
for  spatially  distributed  networked  sensors. 

The  averaging  multisensor  intensity  filter  is  derived  by  Bayesian  methods 
in  [111].  It  is  repeated  here  in  outline.  The  Bayesian  derivation  of  the  single 
sensor  intensity  filter  in  Appendix  D  is  a  good  guide  to  the  overall  structure 
of  most  of  the  argument. 

The  key  is  to  exploit  the  PPP  target  model  on  the  augmented  space  Jv  . 
Following  the  lead  of  Eqn.  (D.5)  in  Appendix  D,  the  only  PPP  realizations 
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with  nonzero  likelihood  have  =  YJ}L  \  mk(f)  microtargets.  The  mj.  PPP 
microtargets  are  paired  with  the  m k  sensor  data  points,  so  the  overall  joint 
likelihood  function  is  the  product  of  the  sensor  data  likelihoods  given  the 
microtarget  assignments.  This  product  is  then  summed  over  all  partitions  of 
the  mjt  microtargets  into  parts  of  size  1), . . . , 

The  sum  over  all  partitions  is  the  Bayes  posterior  pdf  on  the  event  space 
£(»S+).  It  is  a  very  complex  sum,  but  it  has  important  structure.  In  particular, 
the  single  target  marginal  pdfs  are  identical,  that  is,  the  integrals  over  all  but 
one  microtarget  state  are  all  the  same.  After  tedious  algebraic  manipulation, 
the  single  target  marginal  pdf  is  seen  to  be 

Pxsed(*)  =  —  E  Lk{Hk{t)\x-,e)fm-X{x),  xeS\  (6.51) 
mk  t=i 

The  mean  field  approximation  is  now  invoked  as  in  Eqn.  (D.13).  Under  this 
approximation,  /^“sed(x)  =  cp^used(x),  where  the  constant  c  >  0  is  estimated. 
From  (6.17)  and  (6.18),  the  measurement  intensity  is 

A^“sed(z)  =  c  J"  pk(z\x)p^used(x)dx/  (6.52) 

so  the  likelihood  function  of  c  given  the  data  sets  £/,(/')  is 


M 

mk(() 

2^(c;4(1),...,4(M))  =  11 

e-fs+cp^(.r)d,  APused (zk(j;€)) 

e=\ 

i=1 

oc  e~cMcmk . 

Setting  the  derivative  with  respect  to  c  to  zero  and  solving  gives  ML  estimate 
Cml  =  yf  ■  The  multisensor  intensity  filter  is  Further  purely 

technical  details  of  the  Bayesian  derivation  provide  little  additional  insight, 
so  they  are  omitted. 

The  multiplication  of  the  conditional  likelihoods  of  the  sensor  data  hap¬ 
pens  at  the  PPP  event  level,  where  the  correct  associations  of  sensor  data  to 
targets  is  assumed  unknown.  The  result  is  that  the  PPP  parameters  —  the  in¬ 
tensity  functions  —  are  averaged,  not  multiplied.  The  multisensor  intensity 
filter  therefore  cannot  reduce  the  area  of  uncertainty  of  the  extracted  target 
point  estimates.  In  other  words,  the  multisensor  intensity  averaging  filter 
cannot  improve  spatial  resolution.  Intuitively,  the  multisensor  filter  achieves 
variance  reduction  in  the  target  count  by  foregoing  spatial  resolution  of  the 
target  point  estimates. 
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6.5.2  Heterogeneous  Sensor  Coverages 

When  the  probability  of  detection  functions  are  not  identical,  the  multisensor 
intensity  filter  description  is  somewhat  more  involved.  At  each  target  state 
x  the  only  sensors  that  are  averaged  are  those  whose  detection  functions 
are  nonzero  at  x.  This  leads  to  a  "quilt-like"  fused  intensity  that  may  have 
discontinuities  at  the  boundaries  of  sensor  detection  coverages. 

The  Bayesian  derivation  of  (6.47)  outlined  above  assumes  that  all  the 
microtargets  of  the  PPP  realizations  can  be  associated  to  any  of  the  M  sensors. 
If,  however,  any  of  these  microtargets  fall  outside  the  coverage  set  of  a  sensor, 
then  the  assignment  is  not  valid.  The  way  around  the  problem  is  to  partition 
the  target  state  space  appropriately. 

The  total  coverage  set 


C  =  C(0  (6.53) 

contains  points  in  target  state  space  that  are  covered  by  at  least  one  sensor. 
Partition  C  into  disjoint,  nonoverlapping  sets  Bp  that  comprise  points  covered 
by  exactly  p  sensors,  p  =  1,  ...,M.  Now  partition  Bp  into  subsets  Bpri,  ...,BPijp 
that  are  covered  by  different  combinations  of  p  sensors.  To  simplify  notation, 
denote  the  sets  {spy}  by  {J5T„, ),  a>  =  1, 2, . . . ,  O.  The  sets  are  disjoint  and  their 
union  is  all  of  C: 


C  =  JljDJlj  =  0  for  i  +  ]■  (6.54) 

No  smaller  number  of  sets  satisfies  (6.54)  and  also  has  the  property  that  each 
set  in  the  partition  is  covered  by  the  same  subset  of  sensors. 

The  overall  multisensor  intensity  filter  operates  on  the  partition  \J\-a\- 
The  assignment  assumptions  of  the  multisensor  intensity  filter  are  satisfied 
in  each  of  the  sets  Thus,  the  overall  multisensor  filter  is 


rFused 

Jk\k 


(x) 


Wa 


£  4(4(0  \x;€) 


fk\k-l(x)  /  ^  ^  -dr, 


(6.55) 


where  !„,)  are  the  indices  of  the  sensors  that  contribute  to  the  coverage 
of  and  is  the  number  of  sensors  that  do  so. 

The  multisensor  intensity  filter  is  thus  a  kind  of  "patchwork"  with  the 
pieces  being  the  sets  of  the  partition.  The  variance  of  the  multisensor 
filter  is  not  the  same  throughout  C  —  the  more  sensors  contribute  to  the 
coverage  of  a  set  in  the  partition,  the  smaller  the  variance  in  that  set. 
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6.6  Historical  Note 

The  PHD  (Probability  Hypothesis  Density)  filter  ([66],  [67]  and  references 
therein)  was  pioneered  by  Mahler  beginning  about  1994.  Mahler  was  ap¬ 
parently  influenced  in  this  work  [66,  p.  1154]  by  an  intriguing  approach  to 
additive  evidence  accrual  proposed  in  1993  by  Stein  and  Winter  [108].  An 
alternative  and  more  insightful  bin-based  derivation  of  the  PHD  filter  was 
discovered  in  2008  by  Erdinc,  Willett,  and  Bar  Shalom  [29]. 

The  intensity  filter  of  Streit  and  Stone  [117]  is  very  similar  to  the  PHD  fil¬ 
ter,  differing  from  it  primarily  in  its  use  of  an  augmented  target  state  space, 
S+ ,  instead  of  birth-death  processes  to  model  target  initiation  and  termina¬ 
tion.  The  PHD  filter  is  recovered  from  the  intensity  filter  by  modifying  the 
posterior  intensity.  This  paper  is  the  source  of  the  Bayesian  approach  to  the 
intensity  filter  given  in  Section  6.1.  The  other  two  approaches  presented  in 
that  section  follow  the  discussion  given  in  [112].  It  draws  on  the  connec¬ 
tions  to  PET  to  gain  intuitive  insight  into  the  interpretation  of  the  PPP  target 
model.  These  approaches  greatly  simplify  the  mathematical  discussion  sur¬ 
rounding  intensity  filters  since  it  builds  on  work  already  presented  in  earlier 
chapters  for  PET  imaging. 

The  multisensor  intensity  filter  was  first  derived  in  2008  by  Streit  [111] 
using  a  rigorous  Bayesian  methodology,  followed  by  the  same  kind  of  PPP 
approximation  as  is  used  in  the  single  sensor  intensity  filter.  The  general  mul¬ 
tisensor  problem  was  presented  for  both  homogeneous  and  heterogeneous 
sensor  coverage.  The  theoretical  and  practice  importance  of  the  variance  re¬ 
duction  property  of  the  averaging  multisensor  filter  was  also  discussed  in 
the  same  paper. 

Mahler  [66]  reports  a  product  form  for  the  multisensor  PHD  filter.  It  is 
unclear  if  the  product  form  estimates  target  count  correctly.  The  problem 
arises  from  the  need  for  each  of  the  sensor-level  integrals  of  intensity  as  well 
as  the  multisensor  integral  of  intensity  to  estimate  the  target  count.  In  any 
event,  the  multisensor  intensity  filter  and  the  multisensor  PHD  filter  take 
quite  different  forms,  and  therefore  are  different  filters. 

The  MMIF  tracking  filter  is  new.  It  was  developed  by  exploring  connec¬ 
tions  between  intensity  filters  and  the  PMHT  (Probabilistic  Multiple  Hy¬ 
pothesis  Tracking)  filter,  a  Gaussian  mixture  (not  Gaussian  sum)  approach  to 
multitarget  tracking  developed  by  Streit  and  Luginbuhl  [113,  114,  115]  that 
dates  to  1993.  These  connections  reveal  the  PPP  underpinnings  of  PMHT. 
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Bayesian  Derivation  of  Intensity  Filters 


The  multitarget  intensity  filter  is  derived  by  Bayesian  methods  in  this  ap¬ 
pendix.  The  posterior  point  process  is  developed  first,  and  then  the  posterior 
point  process  is  approximated  by  a  PPP.  Finally,  the  last  section  discusses  the 
relationship  between  this  method  and  the  "first  moment"approximation  of 
the  posterior  point  process. 

The  steps  of  the  intensity  filter  are  outlined  in  Fig.  6.1.  The  PPP  interpre¬ 
tations  of  these  steps  are  thinning,  approximating  the  Bayes  update  with  a 
PPP,  and  superposition.  The  PPP  at  time  4  is  first  thinned  by  detection.  The 
two  branches  of  the  thinning  are  the  detected  and  undetected  target  PPPs. 
Both  branches  are  important.  Their  information  updates  are  different.  The 
undetected  target  PPP  is  the  lesser  branch.  Its  information  update  is  a  PPP. 
The  detected  target  branch  is  the  main  branch,  and  its  information  update 
comprises  two  key  steps.  Firstly,  the  Bayes  update  of  the  posterior  point 
process  of  3^  on  B(S+)  given  data  up  to  and  including  time  4  is  obtained. 
The  posterior  is  not  a  PPP,  as  is  seen  below  from  the  form  of  its  pdf  in  (D.10). 
Secondly,  the  posterior  point  process  is  approximated  by  a  PPP,  and  a  low 
computational  complexity  expression  for  the  intensity  of  the  approximating 
PPP  is  obtained.  The  two  branches  of  detection  thinning  are  recombined  by 
superposition  to  obtain  the  intensity  filter  update. 


D.l  Posterior  Point  Process 

The  random  variables  Sw^-i,  and  Tmc-i  are  defined  as  in  Section  C. 

The  state  space  of  i|Jc-i  and  1  is  E(S+),  where  E(S~)  is  a  union  of  sets 
defined  as  in  (2.1).  Similarly,  the  event  space  of  Tq^-i  is  S(T),  not  T. 

The  process  is  assumed  to  be  a  PPP,  so  it  is  parameterized  by 

its  intensity  i|jt_i(s),  s  e  <S+.  A  realization  <4  e  &(S+)  of  Hjt-iiJc-i  is  thinned 
by  a  death  probability  function  c4_i(s),  assumed  known,  and  subsequently 
diffused  to  time  4  via  the  single  target  transition  function  V4-1  (y  I  x) .  Thinning 
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and  diffusing  yields  the  predicted  target  PPP  Euk-i-  Its  intensity  is, 

using  (2.83), 


fk \k-l(x)=  f  xPk-i(x\s)(l-dk_1(s))fk_ljk_1(s)ds.  (D.l) 

JS+ 

The  integral  in  (D.l)  is  defined  as  in  (2.97). 

The  point  process  E uk  is  the  sum  of  detected  and  undetected  target  pro¬ 
cesses,  denoted  by  E®k  and  SjrJ, ,  respectively.  They  are  obtained  from  the 
same  realizations  of  3  m- so  they  would  seem  to  be  highly  correlated. 
However,  the  number  of  points  in  the  realization  is  Poisson  distributed,  so 
they  are  actually  independent.  See  Section  2.9. 

The  undetected  target  process  E^k  is  the  predicted  target  PPP  Ek\k_i 

thinned  by  1  -  Pjr’(s),  where  Pj^(s)  is  the  probability  of  detecting  a  target 

at  s.  Thus  S1:1..  is  a  PPP,  and 

k\k 

f&(x)  =  0--Pk(x))fa-i(x)  (D.2) 

is  its  intensity. 

The  detected  target  process  is  the  predicted  target  PPP  Ek\k_k  that  is 
thinned  by  P^(s)  and  subsequently  updated  by  Bayesian  filtering.  Thinning 
yields  the  predicted  PPP  SjP fc_1,  and 

/£_  iW  =  P?(x)fk  lk-i(x)  (D.3) 

is  its  intensity. 

The  predicted  measurement  process  is  obtained  from  via  the 

pdf  of  a  single  point  measurement  z  e  T  conditioned  on  a  target  located  at 
s  e  »S+.  The  quantity  pk(z  |  (p)  is  the  likelihood  of  z  if  it  is  a  false  alarm.  See 
Section  2.12.  Thus,  Tk \k_i  is  a  PPP  on  T  and 

Afc|fc-l(z)  =  f  p*:(z  | s) Pk  (s)  fk\k-\ (s) ds ,  (D.4) 

Js+ 


is  its  intensity. 

The  measurement  set  is  ujt  =  {//?,  {z\, . . .  ,zm),  where  z(  e  T .  The  conditional 
pdf  of  Ujt  is  defined  for  arbitrary  target  realizations  £*  =  («,  {i'i ,  ...,xn})  e 
£(»S+).  All  the  points  x,  of  whether  they  are  a  true  target  (xj  e  IR";,:)  or 
are  clutter  (x;  =  O),  generate  a  measurement  so  that  only  when  m  =  n  is 
the  measurement  likelihood  non-zero.  The  correct  assignment  of  point  mea¬ 
surements  to  targets  in  i ^  is  unknown.  All  such  assignments  are  equally 
probable,  so  the  pdf  averages  over  all  possible  assignments  of  data  to  false 
alarms  and  targets.  Because  <ft  is  a  target  state,  the  measurement  pdf  is 
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251 


PYk\Bk(Vk\Zk) 


ni\^lo&Sym(tn)  Y\j=i  Pk{Zo(j)\Xj)/  ^71  11 

0,  m  +  n, 


(D.5) 


where  Sym (in)  is  the  set  of  all  permutations  on  the  integers  {1, 2, ,  m). 

The  lower  branch  of  (D.5)  is  a  consequence  of  the  "at  most  one  measure¬ 
ment  per  target" rule  together  with  the  augmented  target  state  space  S+. 
To  elaborate,  the  points  in  a  realization  E,  of  the  detected  target  PPP  are  tar¬ 
gets,  some  of  which  have  state  <p.  The  augmented  state  space  accommodates 
clutter  measurements  by  using  targets  in  (p,  so  only  realizations  with  m  =  n 
points  have  nonzero  probability. 

The  posterior  pdf  of  Sj^.  on  £(»S+)  is,  from  (C.5), 


Pk\k(£,k)  =  PrklSk(Vk\£k) 


Pk\k-l(£k) 

nk\k-l(vk) 


(D.6) 


The  pdf's  pjt|fc_i(4)  and  nm^(vk)  of  Sjj^  and  Tk |fc_i  are  given  in  terms  of 
their  intensity  functions  using  (2.12): 


Pk\k-i(Zk)  =  £ exp £  f°_,(s) dsj  f[  f^k_,(xj)  (D.7) 

nk\k-i(vk)  =  £  exp  (z)  dz  j  Aklk^(z j).  (D.8) 

From  (D.3)  and  (D.4), 

^fkjk-i^^  =  £Am-i(z)dz ■  (D.9) 

Substituting  (D.7),  (D.8),  and  (D.5)  into  (D.6)  and  using  obvious  properties 
of  permutations  gives  the  posterior  pdf  of 


Pklk(Zk)  = 


e  n 

(jeSym(m)  j=  1 


Pk(Za(j )  I  Xj ) fk\k-l(Xj) 

Ak\k-l(za{j)) 


(D.10) 


If  ck  does  not  contain  exactly  in  points,  then  pk\k(Ek)  =  0.  Conditioning  SH, 
on  m  points  gives  the  pdf  of  the  points  of  the  posterior  process  as 


Pk\k{Xl/  ■  ■  ■ ,  Xm) 


1 

ml 


m 


e  n 

creSym(m)  j= 1 


Pk(za(j)  I  Xj)  pD(Xj)  fk\k-l(Xj) 
Ak\k-l(za(j)) 


(D.ll) 


The  pdf  (D.ll)  holds  for  Xj  e  »S+,  j  =  1, ...,  m. 
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D.2  PPP  Approximation 


The  pdf  of  the  posterior  point  process  is  clearly  not  that  of  a  PPP.  This 
causes  a  problem  for  the  recursion.  One  way  around  it  is  to  approximate  SPk 
by  a  PPP  and  recursively  update  the  intensity  of  the  PPP  approximation. 

The  pdf  pk |fc(xi,  =  pk\k(xa(i)'  x<j(»<))  for  a11  a  e  Sym(m);  therefore, 

integrating  it  over  all  of  arguments  except,  say,  the  Pth  argument  gives  the 
same  result  regardless  of  the  choice  of  t.  The  form  of  the  "single  target 
marginal"  is,  using  (D.4), 


Jr*  r*  m 

•••I  Pk\k(xi,  ■■■,xm)Y\dxl 
S+  J  S+ 


1  r  m 

=  i  y  n 

”Z!  creSyin(m)  d(5+)m_1 


Yldx‘ 


(=i 

m 


i  m 

=  s  E  E 


r= 1  (jGSym(m) 
and  o(£)=r 


Pk(za(e)  \x{)P®  (x()  fk\k-i(xi) 

^■k\k-l{za(j)) 


iE 

ra  ' 

r= 1 


Pk(Zr\X{)P%(x{)fkpc-1{x{) 

^■k\k-l{zr) 


(D.12) 


This  identity  holds  for  arbitrary  X[  e  *S+ . 

The  joint  conditional  pdf  is  approximated  by  the  product  of  its  marginal 
pdf's: 


m 

Pk\k(x\, .  ..,xm)  ~  Yl  Pk\k(xj).  (D.13) 

/= i 

The  product  approximation  is  called  a  mean  field  approximation  in  the 
machine  learning  community  [50,  pp.  35-36].  Both  sides  of  (D.13)  integrate 
to  one. 

The  marginal  pdf  is  proportional  to  the  intensity  of  the  approximating 
PPP.  Let  f^k(x)  =  cpk\k(x)  be  the  intensity.  The  likelihood  function  of  the 
unknown  constant  c  is 


1  >» 

J?(c\£k)  =  -e~ycP^ds  H{cpm(xjj)  oc  e~ccm 

h  i 

The  maximum  likelihood  estimate  is  Cml  =  Tit,  so  that 
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m 

= y 


r=  1 


Pk(Zr\x)P°(x)fk\k-l(x) 

Aqfc-l(Zr) 


is  the  intensity  of  the  approximating  PPP. 


253 

(D.14) 


Altogether  Now 

The  PPP  approximation  to  the  point  process  Ek\k  is  the  sum  of  the  undetected 
target  PPP  E^k  and  the  PPP  that  approximates  the  detected  target  process 

E1?., .  Hence, 

k\k  ' 


l-pfw+y; 

r=l 


pk(zr\x)P°(x) 

i(zr) 


fk\k-l(x) 


(D.15) 


is  the  updated  intensity  of  the  PPP  approximation  to  Ek\k. 

The  intensity  filter  comprises  equations  (D.l),  (D.4),  and  (D.15).  The  first 
two  equations  are  more  insightful  when  written  in  traditional  notation.  From 
(D.l), 


fk\k-\  M  =  h(x)  +  J ^  %-i(xls)(l  -  4-i(s))/fc-t|fc-1(s)ds, 
where  the  predicted  target  birth  intensity  is 


h(x)  =  Wk^(x\<f>)(l  ~  4-i(^))A-i|fc-i(^)- 


Also,  from  (D.4), 


A*|jt— i(z) 


Pk(z\s)P°(s)fk\k-i(s)  ds. 


where 


A  k(z)  =  Pk(z\(p)Pk((p)fk\k-i((p) 


(D.16) 


(D.17) 


(D.18) 


(D.19) 


is  the  predicted  measurement  clutter  intensity. 

The  above  derivation  of  the  intensity  and  PHD  filters  was  first  given 
in  [117].  A  more  intuitive  "physical  space"  approach  is  given  by  [29].  An 
analogous  derivation  for  multisensor  multitarget  intensity  filter  is  given  in 
[111]. 
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D.3  First  Moment  Intensity  and  Janossy  Densities 


An  alternative  method  is  often  used  in  the  literature  to  obtain  the  intensity 
function  (D.14)  for  the  detected  target  posterior  point  process.  The  process 
is  not  a  PPP,  but  it  is  a  finite  point  process  because  its  realizations 
contain  exactly  m  points.  The  general  machinery  of  finite  point  processes  is 
thus  applicable  to  An  excellent  reference  for  general  finite  point  process 
theory  is  [16]. 

Let  £  =  (N,  /YIN)  denote  a  realization  of  SR,  where  N  is  the  number  of 
points  and  X\N  is  the  point  set.  From  [16,  Sect.  5.3],  the  Janossy  probability 
density  of  a  finite  point  process  is  defined  by 

jn(x l,---, Xn)  =  PNin)px\N (fri,  •  ■  • , xn) \ n)  for  n  =0,1,2,....  (D.20) 


Janossy  densities  were  encountered  (but  left  named)  early  in  Chapter  2,  Eqn. 
(2.10).  Using  the  standard  argument  list  as  in  (2.13)  gives 

/„(*!,  ...,*„)  =  n\pN(n)px\N(xi,---,xn\n)  for  n  =  0,1,2,... .  (D.21) 

Intuitively,  from  [16,  p.  125], 


jn(x =  Pr 


Exactly  n  points  in  a  realization 
with  one  point  in  each  infinitesimal 
[xj  +  dxj),  i  =  !,...,« 


(D.22) 


Now,  for  the  finite  point  process  ,  pjy(m)  =  1  and  pN(n)  =  0  if  n  i=  m, 
so  only  one  of  the  Janossy  functions  is  nonzero.  The  Janossy  densities  are 


jn  (xlr...,xn) 


m'.pk \k(xi,...,xm), 

0, 


if  n  =  m, 
if  n  +  m, 


(D.23) 


where  pk\k(xi,  •  ■  ■ ,  xm)  is  the  posterior  pdf  given  by  (D.ll).  The  first  moment 
intensity  is  denoted  in  [16]  by  ni  \  (x).  From  [16,  Lemma  5.4.III],  it  is  given  in 
terms  of  the  Janossy  density  functions  by 


mix)  =  >  , 


E  h  L-L^1 . 

n= 0  ‘ 


)  dxj  ■  ■  ■  dx„ 


From  (D.23),  only  the  term  n  =  m  —  1  is  nonzero,  so  that 

mix)  =  (m  ]_  J f  ■  ■  m\pk\kix,  Xi, ..  .,xm_i)  dxi  -  dxm_i . 


(D.24) 


(D.25) 


The  integral  (D.25)  is  exactly  m  times  the  integral  in  Eqn.  (D.12),  so  the  first 
moment  approximation  to  is  identical  to  the  intensity  (D.14). 
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